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ABSTRACT 

We describe a numerical algorithm to integrate the equations of radiation magnetohydrodynamics 
in multidimensions using Godunov methods. This algorithm solves the radiation moment equations 
in the mixed frame, without invoking any diffusion-like approximations. The moment equations are 
closed using a variable Eddington tensor whose components are calculated from a formal solution of 
the transfer equation at a large number of angles using the method of short characteristics. We use a 
comprehensive test suite to verify the algorithm, including convergence tests of radiation-modified lin- 
ear acoustic and magnetosonic waves, the structure of radiation modified shocks, and two-dimensional 
tests of photon bubble instability and the ablation of dense clouds by an intense radiation field. These 
tests cover a very wide range of regimes, including both optically thick and thin flows, and ratios of 
the radiation to gas pressure of at least 10 -4 to 10 4 . Across most of the parameter space, we find 
the method is accurate. However, the tests also reveal there are regimes where the method needs 
improvement, for example when both the radiation pressure and absorption opacity are very large. 
We suggest modifications to the algorithm that will improve accuracy in this case. We discuss the 
advantages of this method over those based on flux-limited diffusion. In particular, we find the method 
is not only substantially more accurate, but often no more expensive than the diffusion approximation 
for our intended applications. 

Subject headings: (magnetohydrodynamics:) MHD — methods: numerical — radiative transfer 



1. INTRODUCTION 

Moving fluids can absorb, emit and scatter photons, and via these processes energy and momentum are exchanged 
between the radiation field and the rest of the flow. When the fluxes of energy and momentum carried by photons 
are significant compared to those carried by particles or magnetic field, the fluid dynamics will be significantly af- 
fected (and perhaps even controlled) by the radiation field. Radiation has been realized to play an important role 
in th e dynamics of many diff e rent astrophysical sys t ems, such as the f ormation of stars in diff erent environments 
/e.g.. iMcKee fc Ostrikerl 120071: IKrumholz et all 120091 : lOffner et al.1 [20091: IjT ang fc Goodman! l20lTh . supernovae (e.g. 



lArnett et al.l 19891: Uanka et al 20071; iNordhaus et al.ll2010j). acc retion flows around supermassive massive black holes 



(e.g.. IShakura fc Sunvaevl 11975 : |Hirose et alll2009t ISpruitll2010l). an d the evolution of galaxies with feedback from a 



central massive black hole (e.g.. ICiotti fc Ostrikerll2007l : lProgall2007l ). to name but a few. 

The dynamics of radiating fluids can be divided into quite different regimes, depending on whether the photons dom- 
inate energy and/or momentum transport, and whether the optical depth across the region of inte rest is small or large 
(an e xcell ent introduct ion to all aspects of radiation hydrodynamics is given in the monographs of Miha las fc Mihalasl 
119841 and ICastorll2004l ). In this work, we are motivated by the properties of accretion flows around compact objects, 
the dynamics of accretion disk boundary layers, and the production of radiation-driven winds and outflows from stars, 
accretion disks, and galaxies. In these systems, radiation often dominates both the energy and momentum transport, 
and the optical depth can vary widely between diffe rent regions in the flow . Since magnetohydrodynamic (MHD) 
effects are crucial in accretion flows (see the review iBalbus fe Hawlevl 119981 and references therein) , they must be 
included. Thus, for these applications the appropriate equations are the those of gas dynamics, Maxwell's equations, 
and the radiation transfer equation. Our goal is to develop numerical algorithms to solve this system of equations in 
order to investigate our motivating applications. 

Although the basic equations of radiation MHD are well known, there is still considerable uncertainty regarding 
even fundamental issues of how to solve them. For example, it is not clear whether to treat the ra diation field in the 
co-moving (fluid) or laboratory frame in order to obtain the simplest and most accurate solutions . ICastorl (12009!) has 
provi ded a detailed critique of both approaches. In this work, we adopt the mixed frame (e.g., I Mihalas fc Mihalasl 
1984), in which both the radiation and fluid variables are integrated in the laboratory frame, but the radiation- material 
interaction terms are treated in the co-moving frame, with 0(v/c) expansions used to transform these terms back to 
the lab frame. A proper accounting of all 0(v/c) terms in these t ransformations is necessary in order to correctly 
account for all effects (e.g.. lLowrie et aL1ll999l : IKrumholz et al.| [2007). Although there are limitations to this approach 
(it is not suitable for treating line transport), we find it is well suited to our applications. 

Rather than integrating the time-dependent transfer equation directly (which, even in the case of frequency inde- 
pendent transport, requires solving an integro-diffcrcntial equation in 6 dimensions), we instead integrate a hierarchy 
of angular moments, using a variable Eddington tensor (VET) that relates the zeroth and second moments in order to 
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close the system of equations. The VET is computed directly from a formal solution of the time-independent transfer 
equati on at every time s t ep. This approach has b een implemented previously in both 2D and 3D yersions of the ZEUS 
code (|Stone et al.lll992t lHaves fc Normanll2003D . the adaptive mesh r efinement code TITAN (|Gehmevr fe Mihalasl 
fl99l as well as a Soften Lagrangian Hydrodynamic code described in iGnedin fc Abel ([2001). However, because of 
the expense of solving the transfer equation in multidimensions, the VET method has not been widely used for astro- 
physical applications. Advances in algorithms and computer hardware now make the VET feasible even in 3D, and 
as we show in this paper, the VET method has considerable advantages over other methods. In a companion paper 
(jDavis et al.ll2012f) . we provide a detailed description of the algorithms we have implemented to solve the transfer 
equation. 

A popular alternative to the VET method for closing the radiat ion moment equations is to adopt the flux-limited 
diffusion (FLD) approximation (e.g.. lLevermore fc Pomranind [T981). In this approach, the radiation flux is assumed to 
be in the direction of the gradient of the radiation energy density, with a value that is limited to prevent superluminal 
transport. In fact, the majority of multidimensional radiation hydrodynamic codes used in astrophysics to date adopt 
FLD fe.g..lTurner fc Stonell200lUKrumholz et al.|[2007t iGittings et al.ll2008tlSwestv fc Mvrall2009l ; Ivan der Hoist et al.l 
I201H ICommercon et al.ll2011t IZhang et al.ll201lD. However, there are a num ber of well-known and potentially serious 
limitations to FLD ('e.g.. iMihalas fc Mihalas 19841 lHaves fc No rman 20031). For example, the fact that the direction 
of the flux is assumed rather than computed from the transfer equation can introduce serious error in optically thin 
regions, such as the inability to follow shadows, or a net force on the fluid which is in the wrong direction. Dropping 
the time-dependent term in the radiation flux eliminates radiation inertia, which can be important when the radiation 
field carries a significant amount of the total momentum in the flow. Finally, in FLD the off-diagonal components 
of the radiation pressure tensor are dropped, which eliminates effects such as radiation viscosity. The VET method 
overcomes all of these limitations. We compare and contrast these two methods throughout this paper. 

A significant difference between the algorithm developed in this paper and previous implementations of the VET 
method is the algorithm used to integrate the hydrodynamic (as well as the MHD) equations. Most previous imple- 
mentations of VET used methods based on op erator splitting, wi th artificial viscosity required for shock capturing as 
implemented in, for example, the ZEUS code ([Stone et al.1 119921 ) . In this work we combine the VET approach with 
a dimensi onally-unsplit hig her-order Godunov method for hydrodynamics and MHD (as implemented in the Athena 
code, see lStone et al.l [2008). By adopting a Riemann solver to compute fluxes of the conserved variables, Godunov 
methods do not require any artificial viscosity for shock capturing. However, a significant challenge to this approach 
is how to treat the stiff source terms associated with the interaction of the radiation and material. With Godunov 
methods, simply splitting these terms from the flux differenc es (the simplest and m ost often adopted approach) is 
known to be problematic when t he source terms are stiff ( e.g., lLowrie fc Mo rel 2001). Instead in this work we adopt 
the modified Godu nov method o f Miniat Tfc Colellal J2007), which provides a stable second-order accurate algorithm. 
In a previous paper ISekora fc Stond (|2010l ) (hereafter SS10), we introduced a one-dimensional version of this algorithm. 
This paper represents the multidimensional extension of the SS10 algorithm. 

In addition to the modified Godunov method, a second crucial ingredient to the algorithm is the method by which 
the VET is computed. We use a formal solution of the radiative transfer (RT) equation based on the method of short 
characteristics, and including multi-frequency, scattering, and non-LTE effects. Angular quadratures of the specific 
intensity are then used to compute the VET from first principles . A complete desc ription of our algorithm for solving 
the RT equation, including tests, is given in a companion paper ([Davis et al.ll2012t ). The RT module can also be used 
to compute images and spectra of MHD simulations for diagnostic purposes, and can even by used to compute the 
radiation source terms in the MHD equations for problems where only energy transport via photons is important. In 
effect, the algorithm described in this paper is a marriage between the modified Godunov method of SS10 (extended 
to multidimensions) and a modern non-LTE stellar atmospheres code to compute the VET. 

Of course, the final measure of any algorithm is provided by quantitative testing. Unfortunately, there are very few 
analytic solutions available in radiation hydrodynamics and MHD useful for such tests. In this paper, we introduce 
a comprehensive test suite that includes error convergence tests of radiation-modified acoustic waves in a variety of 
regimes, quantitative comparison of the structure of radiating shocks in the non-equilibrium diffusion limit to semi- 
analytic solutions over a wide range of Mach numbers, quantitative comparison to previous numerical solutions of the 
structure of subcritical radiating shocks computed with full-transport, the linear growth rate and nonlinear regime of 
the photon bubble instability, and shadowing and ablation of dense clumps by an intense radiation field. It is likely 
that many of these tests will be useful for others developing radiation hydrodynamic codes. 

The outline of this paper is as follows. In the next section, we catalog the equations of motion we solve. In section[3l 
we describe in detail our numerical algorithms, highlighting the extensions we have made to the ID method described 
by SS10. In section [H we discuss the importance of exact energy conservation, and present test problems to measure 
the energy error. We present the results of our comprehensive test suite in section [6] Finally, we compare and contrast 
our numerical algorithms to those based on the diffusion approximation in section [71 and summarize in section [51 

2. EQUATIONS 

Following SS10, we write the equations of radiation MHD in the mixed frame, and solve the ra diation and materia l 
energy and momentum conservation laws se parately, including the Q(v /c) source terms as given bv lLowrie" "eTa!1(fi99"9l) . 
Similar equations have also been derived bv lKrumholz et al.l (120071) in th e flux-limited diffusion approximation, and the 
source terms used here are identical to those in lKrumholz et all (|2007f ) to 0(v/c). Some import ant (v/c) 2 terms are 
also included in our formula, such as the advective flux of radiation enthalpy. See the discussions in lLowrie et al.l (|1999f ) 



3 



on the importance of keeping these terms. We find it convenient to solve the system of equations in a dimensionless 
form by adopting the two ratios C = c/ao and P = a r T^ /Pq. Here ao, To, and Pq are the characteristic values of sound 
speed, gas temperature, and gas pressure respectively, and a r is the radiation constant. Thus, C is the dimensionless 
speed of light, and P the dimensionless radiation pressure. With these scalings, the units of the radiation energy 
density and flux are a r T^ and ca r Tg respectively. 

To further simplify the equations, we assume frequency-independent (gray) opacities, and local thermodynamic 
equilibrium (LTE). The scattering opacity is assumed to be isotropic and coherent in the co- moving frame. Thus, 
we only consider the equations for frequency-integrated quantities, and do not need to distinguish between Planck 
and flux-mean opacities. Extension of our m ethod to frequency dependent transport problems, for example using 
multigroup methods fe.g.. IVavtet et aI1l2011f) . is straightforward but beyond the scope of this paper. The equations 
of radiation MHD are then: 

! + v.(H= 0) 

+ V • (pvv BB + P*) = -PSr(P), 

at 

— + V-[(E + P*)v - B(B ■ v)] = -FCS r (E), 
at 

OB 

— V x (v x B)=0. 

dt 



dE r 

~dT 
dF r 



CV-F r = CS r (E), 

•CV-P r = CS r (P), (1) 



where the source terms are, 



S r (P) - -a t ( F r - vEr V" Pr ) + a a J(T 4 - E r ), 



v ( vE + v ' P \ 

S r (E) =a a (T 4 - Er) + (a a -<*.)■£• \F r r — -J ■ (2) 

In the above, p is density, P* = (P + B 2 /2)\ (with I the unit tensor), a a and a s are the absorption and scattering 
opacities 1 (which can be functions of both density and temperature), and the magnetic permeability jj, = 1. The total 
gas energy density is 

E = E g + l -pv 2 + ^ 1 (3) 

where E g is the internal gas energy density. We adopt an equation of state for an ideal gas with adiabatic index 7, 
thus Eg = P/(7 — 1) for 7 7^ 1 and T = P/RidcaiP, where i?idoai is the ideal gas constant. The radiation pressure P r 
is related to the radiation energy density E r by the closure relation 

P r = fE r . (4) 

where f is the VET. It is straightforward to convert the dimensionless radiation MHD equations given above to their 
more familiar dimensional form by setting P = 1, and replacing C with the speed of light c, T 4 with a r T 4 , and F r 
with F r /c. 

The VET used to close the hierarchy of radiation moment equations is calculated directly from angular quadratures 
of the frequency averaged specific intensity I r 

P r § I r nhdu> 

~E r = §I r du ' (5) 

where uj is solid angle and n is the unit vector. The specific intensity I r is calculated from a formal solution of the 
time-independent transfer equation 

= Kt (S - i r ), (6) 

where n t is the total specific opacity and S the source function. In a companion paper (jDavis et al.ll2012h . we describe 
in detail the algorithm we use to solve the transfer equation [51 which is based on the method of short characteristics. 
In section 13.41 we describe how the angular quadratures of the specific intensity returned by the transfer solver are 
computed to give the Eddington tensor. 

1 Flux mean, Planck mean and energy mean opacities are treated as the same here. But they can be easily extended to be different 
values. 
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In the above equations, radiation quantities are always defined in the Eulerian (lab) frame. To order 0(1/C), the 
co-moving ra diation energy density E r q and radiation flux F r o are related to the lab frame values E r and F r by (e.g., 
lCastorll200l 

v 

E r ,o = E r — 2 — • F r , 

F rj0 = F r --(vE r + vP r ), 

V 

E r = E r fl + 2 — • F rj o, 

F r = F rfi + - {vE rfl + v • P r , ) . (7) 

The source terms (equ ation in t h e mix ed frame representation were originally developed by iMihalas fc Kleinl 
(1982), and extended by lLowrie et al.l Jl999) to include an extra term (a a — a s )v ■ (vE r + v ■ P r )/C 2 in the energy 
source term S r (E). This term is necessary to ensure the correct thermal equilibriu m state in moving fluids, and 
is especially important when scattering opacity is dominant (see a full discussion in lLowrie et al.l |1999( ) . Further 
discussion of the physical interpretation of the source terms can be found in Scction [7.1l 

We do not solve the equations in strong conservation form, which means that total (radiation plus gas) energy 
and momentum are not conserved exactly (to round-off error). However, exact conservation is easily enforced only 
for explicit algorithms, in which case the time step is limited by the light crossing time across a cell. Once implicit 
differencing methods are adopted, then conservation to round-off error is generally not possible due to the much larger 
error tolerance used when inverting the matrix representing the difference equations. Since the use of implicit methods 
is crucial in our application domain, we choose to solve the radiation and material conservation laws with source terms 
as accurately as possible, and monitor the energy error as a diagnostic. The issue of the importance of exact energy 
conservation is discussed in more detail, along with results of tests of energy conservation in our methods, in section 

m 

3. NUMERICAL ALGORITHM 

In a previous paper (SS10) we presented a one-dimensional modified Godunov algorithm for radiation hydrodynamics. 
In this paper our focus will be on the additional extensions and improvements to the SS10 method required in 
multi dimensions. The numerical algorithms described in this paper have been implemented in the Athena MHD 
code (|Stone et al J 12008). Athena implements higher-order reconstruction in the primitive variables, an extension of 
a dimensionally unsplit integrator to MHD, the constrained transport (CT) algorithm to enforce the divergence-free 
constraint on the magnetic field, and a variety of approximate and exact Rieman n solvers. Since comprehensive 
descriptions of the algorithms in Athena has been given previously (|Stone et al.| [2008) , in this section we only describe 
the extensions to these algorithms necessary for radiation MHD. 

One challenge to any numerical algorithm for radiation MHD is the large range of timescales. In our applications, 
we are interested in evolving the system on the sound crossing time, which can be many orders of magnitude larger 
than the light crossing time. Thus, implicit differencing methods are essential. In this work, to improve both stability 
and accuracy, we split an implicit solution of the radiation subsystem from a modified explicit update of the rest of 
the equations. That is, we update the radiation energy density E r and flux F r by solving the moment equations 

dE 

-^+CV -F r = CS r (E), 
OF 

-^+CV-P r = CS r (P). (8) 

using fully implicit, backward Euler differencing. As discussed in SS10, higher-order implicit time integration schemes 
can lead to oscillatory solutions with large time steps. Thus, to ensure a non-oscillatory method, we restrict the 
update to first-order backward Euler. During this step, the gas variables in the source terms S r (E) and S r (P) are 
held constant. 

On the other hand, the gas quantities are updated by solving the ideal MHD equations using a time-explicit modified 
Godunov algorithm for the stiff source terms 

| + v.(H=o, 

+ V • (pvv BB + P*) = -PSp(P), 

at 

f) W 

— + V-{(E + P*)v - B(B ■ v)} = -PCS r (E), 
at 

— -Vx(vxB) = 0. (9) 

In this step, the radiation variables are held constant. Note that the source terms themselves arc not split from the 
flux-divergence terms. This is critical for achieving stable and accurate solutions when the source terms are stiff. 
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The order in which we do these updates is arbitrary. In most cases, we update the gas dynamical variables (that 
is, equations [9]) first. However, when the the radiation pressure is completely negligible, we have found switching the 
order of the update is more robust. 

3.1. Basic Steps in the Algorithm 
To begin with, it is useful to summarize all of the individual steps in the overall algorithm: 

Step 1. — Using the gas variables to compute the opacities and source function, solve the transfer equation over a 
large set of angles using short characteristics, and compute the Eddington tensor f . 

Step 2. — Reconstruct the left- and right-states in the primitive variables at cell interfaces along each of the x, y, 
and z directions independently. We reconstruct the primitive variables, instead of the conserved variables as in SS10, 
since we have found it to be less oscillatory and more accurate. Either second- or third-order spatial reconstruction is 
possible. 

Step 3. — Add the radiation energy and momentum source terms to the left and right states (see equations 79 and 
80 in SS10). Source terms for the left (right) states are calculated by using cell-centered quantities from the cell to 
the left (right) of the interface, and use the modified Godunov method. 

Step 4- — Compute ID fluxes of the conserved variables with the appropriate Riemann solver. For radiation hydro- 
dynamic simulations, we use the HLLC solver, while for radiation MHD simulations we generally use HLLD. We use 
the adiabatic sound speed (instead of the radiation modified sound speed as was used by SS10) to calculate the fluxes. 
This adds extra dissipation and makes the algorithm more robust. We do not calculate fluxes for E r and F r becaus e 
they are not updated in this step; therefore the Riemann solvers are the same as those described in lStone et al.l {2008). 

Step 5. — For radiation MHD calculations, use constrained transport (CT) algorithm, (step 3 of the Athena 3D 
algorithm) to calculate the electric field at cell corners and update the face centered magnetic field. 

Step 6. — Evolve the left- and right-states at each interface with the transverse flux gradients for half time step At/2, 
as required for the unsplit CTU integrator in Athena. 

Step 7. — For radiation MHD calculations, calculate the cell-centered electric field at time f n +V 2 to use as a reference 
state in CT algorithm (step 6 of the Athena 3D algorithm). 

Step 8. — Calculate new fluxes along each direction with the corrected left- and right-states and the appropriate 
Riemann solver. 

Step 9. — Update the area averaged magnetic field at cell faces from time step nton+1 using CT (step 8 of Athena 
3D algorithm). 

Step 10. — Update the gas variables from time n to n + 1 by adding the divergence of the flux gradient and the 
source terms using the modified Godunov method. 

Step 11. — Estimate the value of radiation energy source term added to the gas total energy using the method 
described in Section r3.3.1l Then update the radiation variables E r , F r using fully implicit, backward Euler differencing 
of the radiation moment equations HI This requires solving a large linear system in multidimensions. Gas variables do 
not change during this step. 

Step 12. — Update the time t n+1 — t n + At, calculate the new time step according to the CFL condition with the 
adiabatic sound speed (fast magnetosonic speed) for radiation hydrodynamics (MHD). 

In comparison to the algorithm introduced in SS10, the primary changes we have made to extend the scheme to 
multidimensions are in Step 1 (compute the Eddington tensor in multidimensions), Step 2 (reconstruction in the 
primitive rather than conserved variables), Step 11 (using an estimate of the radiation energy source term, and solving 
the implicit moment equations in multidimension) , and Step 12 (choose a timestep based on the adiabatic sound 
speed). 

3.2. Modified Godunov Method in Multidimensions 

In this subsection, we describe in detail the modified Godunov method to calculate the flux, especially the way to 
add the stiff source terms to the left and right states. 

3.2.1. Reconstruction Step 

To compute the left- and right-states of the vector of conserved variables re quired to calculat e the fluxes via a 
Riemann solver, we reconstruct the primitive variables as originally proposed bv lMiniati &: Colellal (|2007j ). instead of 
the conserved variables (SS10). In hydrodynamics, the primitive and conserved variables are related via 
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pVx 

P Vy 

PVz 

V E J 



(10) 



(7-l)(i£ + «2+«j»)/2 

The radiation source terms for velocity and pressure are 

S r (v) = -¥S r (P)/p (11) 

Sr(P) - (7 - 1)P« • S r (F) - (7 - l)PCS r (E), 

where S T (P) and S r (E) are given by equation [2] These source terms are the same in the case of radiation MHD, even 
though in this case the definition of total gas energy E includes the magnetic energy. Either second- or third-order 
reconstruction schemes can be used (|Stone et al.ll2008ft . 
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The radiation source terms for the primitive variables must be added to the left- and right-state s for one half time 
step. This requires calculating the propagation operator I to project off the unstable mode (e.g., iMiniati fc Colellal 
[20071 SS10), using the gradient of radiation source terms on the plane of primitive variables Vv^S r (W). As discussed in 
SS10, in most cases it is the energy source term S r (E) (or equivalently S r (P)) that defines the stiffness of the problem. 
The momentum source term can be added as a normal body force. However, in the extreme case of radiation pressure 
so large that P > C, the momentum source term may also become stiff. For example, this can happen in the inner 
regions of an accretion d isk around a supermassive black hole, where C ~ 10 2 while P ~ 10 2 (e.g.. lShakura fc Sunvaevi 
ll973HTurner et al.ll2003ft . For this reason, we keep the leading term in the momentum source term, so Vw5 r (W) can 
be written as 



/ 



V w S r (W) = 













qVy 











s 



p 








OP 

Dp 



(12) 



where Sy = dS r {x)/ dy for any quantity x and y. 
For example, in the ^-direction S^/S^ ~ 0{v/C 2 ) and S%JS* ~ 0(v/C 2 ). To order v/C, we can take S£ = 
= = 0, which significantly simplifies the analysis (a significant advantage of using the primitive variables). 
Then the propagation operator is 
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(13) 



\(l - a)P/p a/ 



where we define the following quantities 



p x = [exp(S£ At/2) - 1] /(S£At/2), 

P y = [exp(S^ At/2) - l] m'At/2), 

P z = [exp(S£ Ai/2) - 1] /{Sll At/2), 
a= [exp(S*^At/2) - l] /(S^At/2). 



(14) 



Note that S 1 ^, S v v v y , S v v : - 0(-a t ¥/C) and when a t P/C < 1, then p x , P y , p z sa 1. In this case the propagation 
operator is reduced to the form used in SS10. When atP/C S> 1, then p x , p y and p z ensure that the algorithm is 
still stable in this regime. With this propagation operator, the source terms added to the left- and right-states are 
0.5At x I(At/2)S P (W). 



3.2.2. Fluxes from the Riemann Solver 

Once the left- and right-states are calculated with the appropriate radiation source terms added, the fluxes of the 
conserved variables in each direction can be calculated from any of the Riemann solvers currently implemented in 
Athena (e.g., HLLC for radiation hydro, and HLLD for radiation MHD). Unlike SS10, the characteristic speed uses 
in those solvers is not the radiation modified sound speed (equation 73 of SS10), but the adiabatic sound speed. We 
have found this is necessary to make the multidimensional algorithm stable at a CFL number of 1.0 in 2D, and 0.5 in 
3D (the same stability limits for the extension of the CTU integrator to MHD without radiation). 

In radiation MHD, the cell-centered electromotive force (EMF) must be calculated as a referen ce state in the CT 
algorithm after we get the above fluxes (step 5 of the 2D integrator and step 6 of 3D integrator in lStone et alJ feOOS). 
In this step, momentum needs to be evolved for a half time step and radiation source terms S r (P) need to be added. 

In the multidimensional CTU integrator, the left- and right-states must be corrected with the transverse flux gra- 
dients. Because the radiation subsystem is only first order accurate, transverse gradie nts of the source t erms do not 
need to be added. This is done in the same way as in Athena (equation 86 and 87 of lStone et al.|[2008h . Using the 
corrected left- and right-states, we call the Riemann solver again to calculate the fluxes, which are then used in the 
following predict-correct step. 
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3.2.3. Updating the conserved variables with the predict- correct scheme 

To achieve second-order accur acy, the conserved variab les in this step are updated with a predict-correct scheme, 
similar to the approach taken bv lMiniati fc Colellal f|2007T ) and SS10. The radiation quantities (energy and flux) have 
already been updated by the implicit solution of the radiation subsystem, so they are kept constant during this step. 
To be consistent with our reconstruction step, the stiffness of the momentum source terms in some regions must also 
be taken into account. 

Given the divergence of the fluxes (computed as described above) V • f™+ 1 / 2 an d the radiation source term S r (U n ), 
a predict solution is estimated as 

U = U n + At(\ - AtVuSriU 11 ))- 1 (s v (U n ) - V • F ,l+1/2 ) , (15) 

where Vf/S r (C/ n ) is the gradient of the radiation source term S T (U) with respect to the conserved variables U n as 



V v S r (U) 



The error in this predicted solution, e(Ai), is estimated as 

At 
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(16) 



e(At) = U n + — (s r (U n ) + S r (ufj - AtV ■ F n+1/2 - U. 



Then the correction step is 



fjn+i = u+ q_ AtVuS r {U n ))- l e(At). 



(17) 
(18) 



At the end of this step, all gas and radiation quantities are updated to time step n + 1, and the entire algorithm can 
be repeated in the next cycle. 

3.3. Radiation Subsystem in Multidimensions 

In this subsection, we describe the extensions to SS10 required to integrate the radiation subsystem equations [8] in 
multidimensions. 



3.3.1. Special Treatment of the energy source term 

The radiation energy density E r and flux F r are updated from time step n to n + 1 based on the updated gas 
quantities at time step n + 1, using first-order backward Euler time differencing to make the method stable with a time 
step much larger than the light crossing time. However, in SS10, the source terms on the right hand side of equations 
[5] were calculated using the temperature at time step n. We have found this can introduce significant error in the 
total energy when the gas and radiation are far from thermal equilibrium, and the time step is much bigger than the 
thermalization time. The source of this error is the assumption that the gas temperature is constant during the update 
and thus the energy source terms added to the gas energy density and radiation energy density are not the same. This 
will be worse if there is continuous heating or cooling source terms and the energy error may accumulate. In reality, the 
gas temperature should evolve simultaneously with the radiation energy density until thermal equilibrium is reached. 

To reduce this energy error, we first estimate the change of gas energy density due to the radiation energy source 
term in the modified Godunov step. Let V • F "+ 1 / 2 be the Riemann flux calculated in Section 13.2.21 The updated 
gas total energy after the Godunov step is E n+1 and the gas total energy at the beginning of the step is E n . Then 
the change of gas energy due to the radiation energy source term in this step dS r (E) can be estimated as 

dS r (E) = E n+1 - (E n - dtV ■ F n+1/2 ). (19) 

To conserve total energy, the energy source term added to the radiation subsystem is then S r (E) = 
K [Ca a (f A - E r ) + (a a - a s )v ■ {F r - (vE r + v ■ P r )/C)J - (1 - K)dS r (E)/(dtP), where the estimated temperature 

T 4 = — (1 / (dtCcTa) + l)dSr(E) /¥ + E™. Here 1Z is a parameter we can choose to get the best energy conservation. In 
practice, we find 1Z = 0.05 can give the best results for the tests we have done. 

With this special treatment of the energy source term in the radiation moment equations, we can reduce the energy 
error significantly, especially for states that are initially far from thermal equilibrium in optically thick regions that are 
evolved with a time step much larger than the thermalization time (in practice, such states are extremely rare in any 
dynamical evolution, and usually occur only if set up specifically in the initial conditions). However, we still cannot 
conserve total energy to round-off error. We present tests in section |4] that show in practice, the error in total energy 
conservation is small. 
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3.3.2. The 3D matrix 

The generalization of the fully implicit, backward Euler update of the radiation moment equations to 3D is straight- 
forward. For this step, the vector of conserved quantities at time step n is U™ ■ k = (E r , F r>x , F ryy ,F r ^ z ) . The radiation 
flux F r ^ x , F r y , F r z are advanced to time step n + 1 by solving 

\HLLE 

i,j-l/2,k 

(20) 

while the radiation energy density E r is updated as 

iHLLE 
ij-l/2,k 

(21) 

where the F HLLK are the vector of fluxes for the conserved quantities at each cell interface computed by a HLLE 
Riemann solver (equation 39 in SS10) and dS r (E) is the estimate energy source term fequationlT9"l). Since the backward 
Euler differencing is only first-order in time, first-order spatial reconstruction is used to compute the left- and right- 
states for the HLLE solver. Moreover, since both the HLLE fluxes and the source term S(U™ +1 ) in equation [20l and 
1211 are linear in the unknowns U"^~ fc , then only one matrix solve is required per time step, which can represent a 

significant savings compared to the implici t solution of nonlinear equations that can result from other splittings (e.g., 
iTurner fe~ Stonc 20 Oil iGonzalez et al.ll2007l) . The matrix of coefficients that must be inverted to solve the linear system 
in this step is given in Appendix [X] 
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3.4. Computing the Eddington tensor 

In order to calculate the VET, we use a formal soluti on of time-in d epend ent transfer equation [6] The methods for 
solving the radiative transfer equation are described in iDavis et al.l (|2012f ). and the reader is referred to that work 
for further details. At each time step, and for each grid cell, the specific intensity I r must be integrated along many 
different rays n. The angular discretization and corresponding quadratures are chosen to cover the un it sphere as 
unifo rmly as possible, but still be invariant under 90 degree r otations about the coordinate axes (e.g.. iBruls et aTl 
1999). Along each ray, the method of short characteristics fe.g.. iMihalas et al.lll978H01son fc Kunaszlll987t) is used to 
integrate the transfer equation across the entire mesh. For multidimensional problems, this requires the interpolation 
of intensities, opacities, and emissivities from (only) neighboring grid zones. Simulations with scattering opacity (non- 
LTE problems) are solved iterative l y with accelerated lambda iteration, with methods similar to those described in 
iTrujillo Bueno fc Fabiani Bendichol (11995T). Intensity boundary conditions and parallelization are handled via ghost 
zones, as described in section 3.5 of IDavis et al.l (|2012[ ). 

During the integration along each angle, the zeroth moment (J) and second moment (K) of the specific intensity 
are summed into running quadratures. This saves having to allocate an array to store I r over all angles. For each 
discreet ray hk, there is a vector of direction cosines (fiok, Hik, l^2k) with yuk = hk ■ Xi and quadrature weights Wk- The 
moments are then given by 

n a -l 

k=0 
na-l 

Kij = ^ w khHikHjk, (23) 

fe=0 

where Ik = I r {nk)- Up to a common factor of 4tt/c, J and K are equivalent to radiation energy density and pressure, 
respectively. Since they are computed by the radiation transfer solver, they will (in general) differ from the E r and 
P r of the integrator. As defined in equation (6), the VET is then simply the ratio of these quadratures f = K/J, 
evaluated in each grid cell. 

4. EVALUATING THE IMPORTANCE OF EXACT ENERGY CONSERVATION 

As discussed in section 13. 3[ our algorithm does not conserve energy exactly (to round-off error) , in part because 
we separate the implicit solution of the radiation subsystem from the modified Godunov update of the material 
conservation laws. In fact, even if the strong conservation form of the equations is adopted, it is in general not possible 
to conserve energy to round-off error with implicit differencing if an iterative method is used to solve the resulting 
linear system to some error criterion. Therefore, our philosophy is to monitor energy conservation as a measure of the 
accuracy of the solution, rather than adopting ad-hoc strategies to enforce conservation. 

We have found that one important modification to the original algorithm described in SS10 that improves energy 
conservation is to estimate the energy source term added in the modified Godunov step and then use the same 
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value of energy source term in the radiation subsystem (section I3.3.1[) . A test which demonstrates the importance 
of this improvement is thermal relaxation in a uniform, stationary medium. Consider an infinite uniform gas with 
density p, temperature T, and constant absorption opacity a a . The radiation energy density E r is also uniform 
everywhere and the radiation flux is zero. If the gas and radiation are not in thermal equilibrium, i.e. E r ^ T 4 
in our dimensionless units, then they will evolve towards the equilibrium state E r — T 4 on the thermalization time 
scale ~ P/(FE r a a C) (e.g., Bla es fc Socra tes 2003). The exact solution is set by the condition of energy conservation 
PE r + pR idcal T/(j-l) = ¥E r + pR ideal f/(^-l). We test this evolution using C = 10 4 ,P = l,p= 1,-R idea i = 1,7 = 5/3 
and <7 a = 100, with different initial T and E r , in a domain of size L = 2 and 128 grid points, and periodic boundary 
conditions. In our algorithm the time step is determined by the CFL condition using the adiabatic sound speed, 
which gives about ~ 5 x 10 -3 , much larger than the thermalization time ~ 10 -6 . This problem has also been used 
to test other codes (e.g., SS10, [Zhang et al1l2011l) . but here we choose extreme values of the parameters in order to 
demonstrate potential problems. These values are almost certainly unrealistic in that no dynamical evolution could 
result in these conditions in any cell, but nevertheless it is useful. 

We show the solutions for two different initial conditions: E r = 100, T = 1 (Figure [l| and E r = 1,T = 100 (Figure 
[2]) respectively. Results with our modifications are shown in the left panels of these figures, while the solutions given 
by the method used in SS10 are shown in the right panels. It is clear the our special treatment of the radiation energy 
source term reduces the error in the energy significantly, especially in the case when E r — T 4 is very large. The right 
columns confirm the analysis in section [3. 3. 1[ in that with the original method proposed by SS10, when the time step is 
much larger than thermalization time, E r approaches T in a single step. Unfortunately, this is the wrong solution in 
the case that E r = 1 and T = 100. Because our algorithm adds the same energy source term to the gas and radiation 
energy density, we conserve total energy much better. In this case, the energy error in our algorithm is determined by 
the tolerance level we set in the matrix solver. These tests also demonstrate stability even when the time step is much 
larger than the thermalization time. 

The tests shown above are done with extreme conditions: a very large time step compared to the thermalization 
time and initial conditions far away from thermal equilibrium. The velocity and radiation flux are always zero in 
the tests. As a test of energy conservation in fully dynamical proble ms, we have studied the conservation of total 
energy for magneto- rotational instability (e.g.. iBalbus fc Hawlev|[T998f) in an unstratified shearing-box simulations of 
a radiation dominated black hole accretion disk. Th e initial condition and shearing periodic boundary condition are 
the same as the fiducial model in lTurner et al.l ([20031 ) for location I (see their table 1). The only difference is that the 
ratio between gas pressure and magnetic pressure is 1600 in our simulation. In our units, the initial parameters are 
p = 1, P = 1, E r = 1, B x = 0,B y = 0,B Z = 0.035. The angular frequency is f2 = 1. The dimensionless parameters 
P = 376.26 and C = 4895.61. A resolution of 32 x 128 x 32 grids is used for a box size L x = 1, L y = 4, L z = 1. The 
gas and radiation field are continuously heated up due to the MRI turbulence. The energy source comes from the 
differential rotation of the disk. We calcu late the difference be tween the work done on the simulation box and the total 
energy change according to equation 8 of lHawlev et al.l (|1995l) . which is the energy error. Over 100 period, the energy 
error is only about 0.6% of the final total energy. If measured with respect to the total work done on the simulation 
box, the energy error is only about 1% for 100 periods. 

5. TESTS OF THE RADIATION SUBSYSTEM 

Our implicit Backward Euler scheme to solve the radiation subsystem (equation |5J is only first-order accurate. One 
might be concerned that a first-order integration algorithm may be too diffusive, especially when there are steep spatial 
gradients. In this section, we provide further tests of the numerical solution of the radiation subsystem, and show that 
it can in fact capture sharp features. 

5.1. Marshak Wave Test 

Evolution of a Marshak wave is a ID non-equilibri um diffusion process originally proposed bv iMarshakl (|1958l ). A 
semi-analytic solution is given bv ISu fc Olsoi] (fl99ll . and we compare the results from our code with this solution. 
This test has also been used by many other authors (Gonzal ez et al.l [20071 : iGittings et al.l [20081 : Ivan der Hoist et all 
1201 It IZhang et alJl20lih . It consists of a cold uniform medium with constant absorption opacity a a . A constant 
radiation flux F™ c is applied at the left boundary at t = and allowed to diffuse through the medium. The cold gas 
is heated up by the radiation field, and eventually the gas and radiation field reach thermal equilibrium. This is not 
a dynamical test, so the hydrodynamic algorithm is not used, and only the radiation subsystem is solved, with the 
gas temperature updated according to equation (124) of SS10. All of the other parameters are the same as in SS10. 
In particular, the time step is limited by the light crossing time of each cell in a domain of size L — 0.1 and 512 grid 
points. As shown in Figure [3l our numerical solution matches the semi-analytic solution extremely well. 

5.2. Tophat Test 

The "tophat" or "crooked pipe" test (e.g.. IGentild 120011 ) is designed to study the propagation of a radiation field 
in a low opacity pipe with a complex shape. Surrounding the pipe is a high opacity material. It is a challenging 
problem because it requires following the propagation of radiation through complex shapes bounded by discontinuities 
in opacities. We find this test useful to show that our first-order Backward-Euler differencing with VET can capture 
sharp gradients in the radiation field, as well as follow the propagatio n of radiation along the pipe properly. 

The problem is initialized following the description in iGentild (|2001[) , except that we solve the proble m in Ca r tesian 
instead of cylindrical coordinates. Because of this difference, we cannot reproduce the solution given bv IGentild ([200l 
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Fig. 1. — Evolution of the radiation and gas energies for relaxation to thermal equilibrium using a ti mestep ~ 10 3 larger than the 
thermalization time. The left panel is the result with our special treatment of the energy source term fsection l3.3,ll l while the right column 
is the result with the original method used by SS10. The red (black) lines are for E r (T 4 ) while the blue lines are the thermal equilibrium 
state. 



exactly. Nonetheless, the test is still useful to demonstrate we achieve qualitatively the same solution, and that the 
VET method can solve this problem properly. 

The size of the simulation domain is (0, 7) x (—2, 2) with a resolution of 512 x 256 grid points. Dense, opaque 
material with density 10 g/cm 3 and opacity <r a = 2000 cm" 1 is located in the following regions: (3,4) x (—1,1), 
(0,2.5) x (-2,-0.5), (0,2.5) x (0.5,2), (4.5,7) x (-2,-0.5), (4.5,7) x (0.5,2), (2.5,4.5) x (-2,-1.5) and (2.5,4.5) x 
(1.5,2.5). The pipe, with density 0.01 g/cm 3 and opacity a a = 0.2 cm -1 , occupies all other regions. The structure 
of the pipe is shown by the black line in the top panel of Figure [4] The dimensionless speed of light is 300. Initially, 
the material has a temperature 0.05 keV everywhere, and the radiation and material temperature are in equilibrium. 
A heating source with a fixed temperature 0.5 keV is located on the left boundary for —0.5 < y < 0.5. All other 
boundary conditions are outflow. We use the short-characteristic module to calculate the VET. An isotropic incoming 
specific intensity is applied only along the left boundary in the region covering the heating source. The incoming 
specific intensity is zero for all other boundaries. During the evolution, velocities are always kept zero so the material 
density does not change. Thus, the modified Godunov step to update the material quantities is not needed in this 
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Fig. 2. — Same as figure [T] but with a different initial condition. Although the initial difference of T 4 — E r is ~ 10 s , our improved 
algorithm relaxes to the equilibrium state accurately. 



test, and we simply evolve the material temperature through the following two equations 

c vP (T n+1 - T n ) = -dtCa a (a r T n+1A - E? +1 ), 

E? +1 - E> r l =dtCa a (a r T n+1A - E? +1 ), (24) 

where the heat capacity c v — 10 15 erg g _1 keV" 1 . Note that we only update material temperature in this step; the 
radiation energy density is unchanged. The radiation energy source term added in this step is dS r (E) = c v p(T n+1 —T n ). 
The radiation energy density and flux are then evolved using our first-order Backward Euler update. We call this 
algorithm method I. 

First, we fix the time step to be 4 x 10 -5 to resolve the light crossing time. Snapshots of the temperature distribution 
at two different times are shown in Figure U] Note that the shadows around the corners of the pipe are captured 
properly. Five probes are placed at (x = 0.25, y — 0), (x — 2.75, y = 0), (x — 3.5, y — 1.25), (x = 4.25, y = 0) and 
(x = 6.75, y = 0) to monitor the change of the temperature in the pipe. The time evol ution of the te mperatures at 
the five points are shown as solid black lines in Figure [SJ Compared with Figure 9 of iGentild (|2001[ ). which shows 
the result calculated with a Monte Carlo method, our first-order backward Euler method with VET gives very similar 
results when the time step is small enough to resolve the light crossing time. For the fifth p oint, no t e that it cools off 
slightly before being heated by the radiation wave, in agreement with the behavior noted bv lGentilel <|2001l ) . The most 




Fig. 3. — Marshak wave test. The filled circles and squares are the semi-analytic solution of Su & Olson (1996), while the solid and 
dashed lines are numerical solutions, for the fourth power of the gas temperature and the radiation energy density respectively. Each pair 
of lines with the same color are at a different time, with time increasing from the bottom to the top pair. 



obvious difference in our solution compared to iGentilej (|2001l) is that the temperature of the first two points increases 
too quickly. This is likely due to the fact that we neglect the light travel time in our short-characteristic module when 
we calculate the VET. Thus the first two points are affected by the hea t ing so urce too early. We also find the evolution 
of the last three points to differ in detail from that shown by IGentilej (|2001| ) , however for these points the difference 
between propagation in cylindrical versus planar (Cartesian) geometry may be important. In our solution, there is 
no geometrical dilution of the flux as it propagates radially, and therefore the heating rate per unit surface area on 
the walls of the pipe is increased. This can affect the detailed temporal evolution at the third through fifth probes. 
Most importantly, note that from Figure 2] our first order implicit update can maintain very sharp gradients in the 
temperature at the walls of the pipe. 

As a second test of our method, we repeat the problem using a time step of 10~ 3 initially, and increasing it by a 
factor of 1.1 at each step until it reaches 1. The result is shown as the dashed line in Figure [S] With this larger time 
step, the error in the temporal evolution of the probes is increased, although the basic evolution is still correct. This 
is consistent with our expectations. If we are not interested following evolution on a light crossing or thermalization 
time, we can use a large time step and still maintain a stable solution with the correct asymptotic behavior. Instead, 
if we want to explore the evolution on these time scales, we must use a small enough time step to resolve them. 

A second method to evolve the material temperature with our algorithms is to add the radiation energy source term 
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Fig. 4. — Temperature of the pipe in the Tophat test at time 0.8 (top) and 9.4 (bottom) respectively. Temperature unit is 0.5 keV. The 
five black dots in each plot are the positions where the temperature probes are placed. The heating source with fixed temperature 0.5 keV 
is located at the left hand side of the pipe. The black line in the top panel is the boundary of the pipe. 



based on the solution of the time- independent transfer equation used to calculate th e VET. Deta i ls of h ow to compute 
the material heating rate from the solution of the transfer equation are given in iDavis et al.l ([2012). We call this 
algorithm method II. The solution computed using this method for the tophat test is shown as red lines in Figure [5j 
The time step is fixed to be 10~ 3 . The result does not change substantiall y when we dec rease the time step further. 
This method also gets a very similar solution as that shown in Figure 9 of lGentild l|200l . Once again, the material 
is heated up too quickly everywhere but the location of the first probe. Again, this is likely due to the fact that the 
short characteristic module solves the time-independent radiation transfer equation, so that we ignore the delay due 
to the finite speed of propagation of radiation down the pipe. The last three points also differ because of our use of a 
planar geometry. Nevertheless, method II still captures the basic evolution correctly. 

Other stringent te st of our radiation transport module used to calculate the VET (such as the beam test), are 
described in detail in lDavis et al.1 (|2012| ). and will not be repeated here. 

6. TESTS OF THE FULL ALGORITHM 

In SS10, tests of various aspects of the one dimensional version of the algorithm used here were presented, and we 
have verified that our extension to multidimensions is able to pass all of these same tests. However, more rigorous 
testing requires solutions to the full system of equations of radiation hydrodynamics in multidimensions. Unfortunately, 
unlike the case of hydrodynamics or MHD, there are very few such solutions available. In the following subsections, 
we describe the results of the suite of test problems we have found most useful. 

6.1. Linear Wave Convergence 

The dispersion relation for linear wave solutions to the radiation hydrodynamic equations in a uniform homogeneous 
background medium, and as suming the Eddington approximation (fix ed E ddington factor f = 1/3 1), have been analyzed 
by a variety of authors (e.g.. lMihalas fc Mihalall984lBogdan et al.|[l99l Uohnson fc Kleinll2010h . T hese authors solve 
the bo undary value problem in which a complex wave number k is found for each real frequency to. iJohnson fc Kleinl 
(2010) have used these solutions to test the Lagrangian radiation hydro dynamic code KULL by driving the waves with 
a time-dependent boundary condition; similar tests were used earlier bv lTurner fc Stone] (|2001l) to test the flux-limited 
diffusion module in the ZEUS code. 
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Fig. 5. — Evolution of temperature for the five probes in the crooked pipe. From left to right, each line corresponds to one of the 
temperature probes labeled in Figure [4] The solid line is the solution from method I with fixed time step 4 X 10 -5 to resolve the light 
crossing time. The dashed line is also the solution from method I with time step 10 -3 to begin with, and increased by a factor of 1.1 at 
each step until the time step reaches 1. The red line is the solution from method II with a fixed time step 10 — 3 . 



In contrast. lLowrie et al.l ([19991) give the dispersion relation for linear waves (again in the Eddington approximation) 
as an initial value problem, that is solving for the imaginary frequency u at fixed real k. Although their analysis 
focused on the properties of linear waves in a moving fluid (they point out that in radiation hydrodynamics, linear 
wave solutions are not Galilean invariant), their approach is ideally suited for code tests. It allows the propagation of 
linear waves to be followed in a periodic domain, free from the complexity of the implementation of driving boundary 
conditions in an Eulerian code, which may in and of itself introduce s purious error into the sol ution s. The dispersion 
relatio n for linear waves in radiation MHD has also been considered by iBlaes &: Socratesl (|2001[ ) and iBlaes fc Socratesl 
(|2003h . We use both the hydrodynamic and MHD linear wave solutions for convergence testing in ID, 2D, and 3D 
in the following subsections. We emphasize that in order to keep the solutions to the dispersion relation analytic so 
that they can be used in the test problems, the Eddington approximation must be adopted. There may be regions 
of parameter space where the Eddington approximation is not appropriate, and the properties of linear waves in this 
regime should be computed using a variable Eddington factor computed self-consistently with flow. However such 
solutions are beyond the scope of this paper. 

It is not possible to simply initialize this test problem using previous solutions for linear waves given in the literature, 
because there may be differences in the frame of reference, and in the order of the source terms kept, in the fundamental 
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system of dynamical equations (our equation [T]) . For this reason, we have rederived the dispersion relation for linear 
waves in the system of equations solved by our code. Even with the Eddington approximation, these solutions are 
non-trivial, requiring finding the roots of a fourth-order polynomial for complex k in hydrodynamics, and a tenth-order 
polynomial for k in MHD. Appendices B and C give the details of our solutions in the case of hydrodynamics and 
MHD respectively, while in Tables 1 through 3 we give the complete eigenvector for linear waves for several parameter 
values for use by others. 

6.1.1. Linear Waves in Radiation Hydrodynamics 

We begin by studying the propagation of linear waves in one dimension with hydrodynamics. We consider a uniform, 
homogeneous medium with background state p = T = P = 1, adiabatic index 7 = 5/3, and zero velocity initially. 
The radiation flux is zero, and the radiation energy density E r = 1 initially for all of the calculations. We assume a 
constant absorption opacity o~ a (perturbations of the opacity will always be second order, which can be neglected), and 
zero scattering opacity. In order to adopt the Eddington approximation, we fix the Eddington tensor to be diagonal 
with components fa — 1/3. The dimensionless speed of light C — 10 4 . We use a domain of size L — 1 with 512 grid 
zones, and periodic boundary conditions for all variables at each edge. We initialize a wave solution (Appendix [B]) by 
adding an eigenmode (see Table 1) with a wavelength equal to the size of the domain i, and an amplitude of 10~ 6 . 

We perform a series of calculations in which we study the effect of varying the radiation to gas pressure ratio by 
varying the dimensionless pressure P, while keeping the initial conditions fixed as above. We also study the effect of 
varying the optical depth per wavelength r a = o~ a L by varying the absorption opacity a a . For each of these calculations, 
we measure the phase velocity of the waves by determining how long it takes a fixed feature in the wave (the density 
maximum) to propagate a distance equal to one wavelength (L = 1). We also measure the damping rate of the wave 
by determining the best fit amplitude to the entire wave profile after one period. We then compare these results to 
linear theory. 

The properties of radiation modified linear acoustic waves vary dramatically with optical depth and pressure ratio. 
Two dimensionless numbers can be used to define various regimes. The first, PCr a , measure the importance of energy 
exchange between the radiation field and the material, while the second, Pr a , measures the importance of momentum 
exchange. When both PCr a <C 1 and Pr a <C 1, energy exchange between the radiation field and material is very 
weak, the momentum of the radiation field can be neglected, so the waves are weakly damped and propagate at nearly 
the adiabatic sound speed. The damping rate increases with optical depth until r a ~ 1; thereafter the radiation 
and material energy densities are strongly coupled, and the damping rate decreases with increasing optical depth. 
When PCr a ^> 1 and P < 1, the radiation and material energy densities are strongly coupled but the radiation 
momentum is still negligible, therefore in this case the damping rate decreases with optical depth until r a ~ 1, and 
increases thereafter. Finally, if PCr a ^ 1 and P > 1, the momentum in the photons becomes important, and the 
damping rate reaches a minimum when Pr a ~ 1 and increase afterwards. To ensure our algorithms are accurate 
over a wide range of regimes, it is important to perform tests that span these dimensionless parameter values. These 
dimensionless quantities a lso clarify a potential shortcoming with numerical algorithm based on the reduced speed of 
light approximation (e.g., iGnedin fc Abelll2~00lD . Any arbitrary reduction in the speed of light will reduce C and the 
above dimensionless numbers correspondingly, which will alter the damping rate and phase velocity of the radiation 
modified acoustic wave. 

In addition to P and C, the Boltzmann number Bo is another dimensionless number w hich is useful to quantify the 
relative importance of radiative and material energy transport in a r adiating flow (e.g.. iMihalas fc Mihal"aslll984D ; it 
is defined to be the ratio between the material enthalpy flux and radiative flux from a free surface. It is related to our 
dimensionless numbers P and C as 

Bo = ¥c^Ty (25) 

where v = v/a$, is the typical flow velocity in units of a$. In the linear wave tests, the typical flow velocity v can 
be replaced with the adiabatic sound speed. When Bo is small, energy transport is dominated by radiation field and 
material energy transport is dominant when Bo is large. 

Figure [5] compares the numerically measured phase velocity and damping rate (shown as stars) for two different 
radiation pressures, and over a wide range of optical depths r a , in comparison to the solution of the linear dispersion 
relation (equation IF31|) shown as solid lines. The parameters are chosen such that for the runs with P = 10~ 4 , the 
dimensionless number PCr a spans 10~ 2 < PCr a < 10 2 , while for the runs with P — 1 both of these limits are a factor 
of 10 4 larger. Based on the discussion above, for P = 10~ 4 we expect the damping rate to increase until r a ~ 1 
and decrease thereafter, while for P = 1 we expect the opposite. The solid lines from the analytic solution of the 
dispersion relation clearly show these trends. In addition, the numerically measured phase velocity agrees with the 
analytic results over the entire range of optical depths for both values of the radiation pressure, while the damping 
rate is accurate in all cases except when P = 1 and r > 1. However, for these parameter values, the damping rate 
is small and the numerical measured damping rates are not converged at the resolution used for this plot. We have 
confirmed that if we increase the resolution, the numerical measured damping rate converges to the analytic values. 
At the same time, we have also found when r a > 100, our algorithm no longer reproduces the correct damping rate. 
We discuss this further, and suggest a solution, below. When P = 10~ 4 , Boltzmann number Bo = 10 and it is 10~ 3 
when P = 1. In the case P = 100, Bo is only to 10~ 5 . 
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Fig. 6. — Phase velocity scaled to the isothermal sound speed (top panel) and damping rate (lower panel) of linear acoustic waves versus 
the optical depth per wavelength for two different radiation pressures. T he st ars represent quantaties measured from simulations, while 
the lines come from the solution to the linear dispersion relation (equation IB II) . All simulations are performed in ID with 512 grid points 
per wavelength. 



In Figure [7J we test linear wave propagation over a wide range of radiation pressures. Once again, the left panel 
compares the numerically measured phase velocity and damping rate (shown as stars) for three different radiation 
pressures, and over a wide range of optical depths t q , in comparison to the solution of the linear dispersion relation 
(equation IBip shown as solid lines. Even though the properties of radiation modified linear acoustic waves vary 
dramatically over this range of optical depths and pressure ratios, there is good agreement between the numerical 
and analytic solutions, except when P = 1 and r a > 1. For these parameters, we also find at higher resolution, our 
algorithm converges to the analytic damping rate, suggesting that the discrepancy is dominated by numerical diffusion 
at this resolution. 

Note that in Figure [7] we do not plot the numerical solution for r = 100 and P = 100 or 0.01. We have found 
the damping rate of linear waves stops converging for these parameter values. Our tests reveal that the problem is 
associated with our implicit solution of the radiation subsystem at a timestep determined by the adiabatic sound speed 
in the material. In this regime, an accurate solution require a timestep which resolves the light crossing time across 
a cell. Note we find this is only an issue for very large absorption opacity: large scattering opacities do not couple 
the radiation and material energy densities. When the optical depth per wavelength due to the absorption opacity 
becomes so large (Pr a ~ C) that the radiation field and the material can be treated like a single fluid, we have found 
it is more accurate to solve a single system of conservations laws for the total (radiation plus material) energy and 
momentum, rather than splitting the radiation subsystem from the material conservation laws. 

As our next test, we measure the convergence rate of linear waves in 3D, for two different values of the problem 
parameters. We use a computational domain of size 21x1x1, with the same number of grid cells per unit length L 
in each direction. We initialize the one dimensional solution inclined to the grid using the same technique to compute 
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Fig. 7. — Left. Phase velocity scaled to the isothermal sound speed (top panel) and damping rate (lower panel) of linear acoustic waves 
versus the optical depth per wavelength for three different radiation pressures. The stars represent quantities measured from simulations, 
while the lines come from the solution to the linear dispersion relation (appendix [B} . All simulations are performed in ID with 512 grid 
points per wavelength. Numerically measured values of both quantaties agrees with linear theory when the absorption optical depth per 
wavelength is smaller than ~ 100. Right. Convergence of LI error with resolution for two different sets of parameters. Both simulations 
are in full 3D. The radiation dominated case converges at first order, while the gas pressure dominated case converges at second order as 
expected (see text). 



the appropriate volume averages of the solution as described in ([Stone et al.ll2008l ). We propagate the wave for one 
period, and then measure the LI error in the solution from 

i 

where is the initial solution, and the sum is taken over all grid points. We repeat this calculation for a variety of 
different numerical resolutions (grid cells per unit length L), and plot the change in the LI error with resolution. The 
result is shown in the right panel of Figure [7] for a radiation dominated, optically thick fluid (P = 10 2 , r a = 10), denoted 
by the solid circles, and for a gas pressure dominated, optically thin fluid (P = 10~ 2 , r a = 10~ 2 ), denoted by the open 
squares. In the radiation dominated case the errors converge close to first order, while in the gas pressure dominated 
case, the errors converge close to second order. This behavior is expected. The errors in the former are dominated 
by the solution to the radiation subsystem, which uses a first-order accurate backward Euler step for stability. Thus, 
we expect the overall rate of convergence should be first order. On the other hand, in the gas pressure dominated 
case, the errors are dominated by the solution to the material conservation laws, which uses a second-order accurate 
modified Godunov step. In this case, the overall rate of convergence should be second order. Of course, it might be 
possible to increase the rate of convergence in the implicit solution of the radiation moment equations by adopting 
higher order implicit differencing. However, limiting such methods to enforce monotonicity can be problematic (e.g., 
SS10). We prefer to adopt unconditionally stable first-order methods for the implicit solver. 

Linear wave solutions wh en only radi a tion e nergy source terms are added with radiation momentum source terms 
neglected are also given by iDavis et al.l (|2012t ) and used to test our radiation transfer module (see Figure 9 of that 
paper) . Because the radiation energy source term is added explicitly, the radiation diffusion mode needs to be resolved 
to make th e code stable, which can limit the time step significantly when Bo < 1. Therefore the algorithm of 
IDavis et al.l ((2012) is most suitable for the regime Bo > 0.01. Because we include all the radiation mom entum and 
energy source terms, the dispersion relations given in this paper differ from Figure 8 of IDavis et al.l (|2012t ). However, 
in regimes of parameter spaces which overlap, we do get very similar results, for example when P = 10~ 4 and Bo = 10 
shown in Figure [6] 

6.1.2. Linear Waves in radiation MHD 

Next, in order to test the accuracy of our MHD algorithms with radiation, we consider the propagation of linear 
modes of radiation modified magnetosonic waves. We do not consider Alfven waves in this subsection since they are 
incompressible and are affected less by radiation than magnetosonic modes. In the case of radiation modified MHD 
waves, the transverse components of the velocity must be included even for ID solutions. However, as implemented our 
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code does not include the transverse velocities in ID during the implicit solution of the radiation moment equations. 
Thus, all of the MHD tests presented here have been performed in 2D using a grid of either 512 x 64, or in some cases 
1024 x 64 cells. 

For these MHD tests, we use a uniform, homogeneous medium with the hydrodynamic and radiation variables set 
to the identical values as were used for the hydrodynamic test (see the previous subsection). The strength of the 
magnetic field is Bq = v/5/3 (which gives a ratio of the Alfven to sound speed of 2). The direction of the magnetic 
field is 45 degrees to the x— axis, and it is confined to the x — y plane (so B z — if). We use a 2D domain of size L x L 
with L — 1 with periodic boundary conditions for all variables at each edge. As before, we initialize a wave solution 
(Appendix [C]) by adding an eigenmode (see Tables 2 or 3) propagating in the x— direction with a wavelength equal 
to then size of the domain L, and an amplitude of 10 -6 . Once again, we perform a series of calculations in which we 
study the effect of varying the radiation to gas pressure ratio by varying the dimensionless pressure P, and the effect 
of varying the optical depth per wavelength r a — a a L by varying the absorption opacity a a . 

The results of our tests for slow magnetosonic waves are shown in Figure while the results for fast magnetosonic 
wave are shown in Figure [9l The left panel in each figure compares the numerically measured values of the phase 
velocity and damping rate with linear theory. Note that the solution to the dispersion relation as a function of 
optical depth and P for both magnetosonic modes is very similar to that for radiation modified acoustic waves in 
hydrodynamics (see Figure [9]). Once again, we find good agreement between our numerical values for the phase 
velocity versus linear theory over this parameter regime. The damping rate also shows good agreement, except when 
P = 100, r a < 1 and P = l,r Q > 0.1. In these cases, the damping rate is so small that our numerical measured 
values are not converged at this resolution. Nevertheless, as Figure [9] demonstrates, our algorithm converges to the 
correct solution for P = 1,t q = 1, albeit at first order. We have confirmed our algorithm also converges for the other 
parameters shown in this figure as well. 

The right panels of Figures [5] and [S] show the convergence rate of the LI error with numerical resolution for the slow 
and fast magnetosonic waves respectively in a fully 3D domain of size 2L x L x L. As before, we measure first-order 
convergence in a parameter regime where the solution to the radiation moment equations dominates the error, and 
second-order convergence in the opposite limit. 




Fig. 8. — Phase velocity scaled to the isothermal sound speed (top panel) and damping rate (lower panel) of linear slow magnetosonic waves 
versus the optical depth per wavelength for three different radiation pressures. The stars represent quantaties measured from simulations, 
while the lines come from the solution to the linear dispersion relation (appendix [Cj. All simulations use an effective resolutions of 512 
grid points per wavelength (except for the two runs with P = 100 and r a < 1, which use 1024 grid points per wavelength). Numerically 
measured values of both quantaties agree with linear theory over the entire range of parameter space, except for very low damping rates, 
where numerical diffusion dominates. Right. Convergence of LI error with resolution for two different sets of parameters. Both simulations 
are in full 3D. The radiation dominated case converges at first order, while the gas pressure dominated case converges at second order, as 
in the case of hydrodynmaics. 



6.2. Radiative shocks in the non- equilibrium diffusion limit 

Shock tubes have long been used as a test of hydrodynamic and MHD codes in the nonlinear regime, since the 
structure of a planar shock can be computed analytically in these cases. However, in radiation hydrodynamics and 
MHD, shock structure is much more complicated, and is difficult to compute analytically. At low Mach numbers A4, 
radiation can diffuse upstream of the shock, heating the gas and forming a smooth precursor in the temperature. As 
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M. increase s, the gas temperature near the shock front can begin to exceed the downstream value, forming a Zel'dovich 
spike (e.g... iZel'Dovich fc Raizcr 1967; Mihalas & Mihalas 1984). This spike is followed by a relaxation region where 
the gas temperature cools to its far downstream value. When a Zel'dovich spike is formed, if the downstream gas 
temperature (after the relaxation region) is larger than the value immediately upstream of the shock (at the end of 
the precursor region), the shock is called subcritical. On the other hand, if the downstream gas temperature after the 
relaxation region is the same as the upstream gas temperature, the shock is termed supercritical. Subcritical shocks are 
formed at lower M than supercritical shocks. Finally, as M. is increased further, the downstream radiation pressure 
can exceed the gas pressure, the Zel'dovich spike disappears, and the solution becomes everywhere smooth again. 

Although many authors have presented numerical solutions to radiating shocks as a t est of their algorithms, the 
lack of analytic solutions inhibits quantitive testing. Recently lLowrie fc Edwards! (|2008f ) have studied the structure 
of radiation modified shocks in the non-equilibrium diffusion limit as a function of Mach number. They give a very 
clear explanation of how these structures can be computed by combining smooth solutions to ordinary differential 
equations with discontinuous jumps determined by the Rankine-Hugoniot relations when needed. These solutions 
not only provide interesting insight into the structure of radiating shocks, but also provide an quantative test of 
time-dependent radiation hydrodynamic codes. 

To test our algorithms, we follow the semi-analytic method described in lLowrie &: Edwards! (|2008l ) (hereafter LE) to 
calculate the shock structure at a variety of Mach numbers, using a non-dimensional pre-shock solution of po = To = 
E r fi — 1 and vq — M. in the upstream state. We then initialize this solution on a ID grid by volume averaging each 
conserved variable to our numerical mesh. The LE solutions are given in the co-moving frame, thus we transform the 
co-moving radiation flux and energy density to the Eulerian frame, using equation [7J In order to make our solutions 
match those presented in LE, we use a dimensionless speed of light C = 1.732 x 10 3 and pressure P = 10~ 4 and 
absorption and scattering opacities of a a = a t — 577.4. The gas temperature T is calculated as T = P/(0.6p) and 
7 = 5/3. The Eddington approximation is used so P r — E r /3. We use a computational domain of size L, where L is 
large enough to capture the upstream precursor and downstream radiative relaxation regions (the size of these regions 
are given by the LE solution itself), with a grid of 1024 cells in all cases (except Ai — 50 where we we use 2048 cells). 
Input boundary conditions at the upstream values are used on the left, and outflow boundary conditions are used 
on the right, with the shock propagating in the negative x-direction. We then let the code evolve this solution for 
several flow crossing times, tf = L/Ai. Ideally, the solution should remain stationary on the mesh. A small shift in 
the position of the shock front can be expected due to truncation error in the averaging of the initial solution to the 
hydrodynamic grid. A serious error in the algorithm or its implementation would be revealed if the code cannot hold 
the input solution. 

We begin by presenting our numerical solution for Ai = 1.05 in Figure ITOl In this case there is no discontinuous 
jump in any variable, but rather only a smooth precursor. Next, Figure fTTI shows the solution for Ai — 3. Now a 
discontinous shock jump has appeared, along with a Zel'dovich spike. For these parameters, this shock is subcritical. 
The inset panel in the plot of gas temperature shows a blow-up of the region near the spike. The numerical solution 
shows agreement with the semi-analytic solution of LE. Next, Figure fT2l shows the solution for Ai = 5. Now there is a 
discontinuous jump in some variables (such as density and velocity), but the gas temperature is continuous apart from 
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the Zel'dovich spike. This solution corresponds to a supercritical shock. Again, the blow-up of the gas temperature 
near the spike region reveals agreement between the numerical and semi-analytic solutions. The position of the spike 
has moved a few cell widths, but has quickly relaxed to a steady structure at the correct shock speed. Finally, Figure 
1 131 shows the solution for M — 50. In this case, the radiation pressure is about 10 times gas pressure in the downstream 
flow. Now all variables are smooth, with no jumps. 
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Fig. 10. — Strucure of a radiation modified shock for Mach number JV[ = 1.05. Blue lines are initial conditions, while black dots are 
numerical results after a few flow crossing times. 



LE have shown that the gas temperature profile near the shock front is very sensitive to small changes in A4 when M. 
is small. To explore whether our algorithms can capture these subtle changes, Figure [14] plots the gas temperature for 
six Mach numbers between A4 = 1.05 and 5. The values are chosen to correspond to figures 4 through 11 of LE, where 
these changes are noted and discussed. This region covers the formation of a Zel'dovich spike, and the increase in the 
amplitude of the spike with increasing M. . It is clear from the figure, in which the numerical solution is compared to 
the corresponding solution from LE, that our algorithm captures each phase accurately. 

Figure [TJ] shows that our algorithm captures the transition from subcritical to supercritical shocks accurately at 
low M., while figure [T3l shows the solution with A4 = 50 (in the stongly radiation pressure dominated regime) is also 
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Fig. 11. — Same as figure [TO] but for M = 3. At this Mach number, a subcritical shock is formed. Inset in the gas temperature panel 
shows details of the Zel'dovich spike. 



captured accurately. However, we have found solutions at M. f» 30 are inaccurate at the resolutions we use. As the 
Mach number increases beyond M. = 5, the width if the Zel'dovich spike compared to the width of the precursor region 
gets smaller, until the spike disappears in the radiation pressure dominated case near M. — 50. We find that unless 
the grid resolution is smaller than the thickness of the Zel'dovich spike, our algorithm will not hold a steady solution. 
For example, at M. = 30 this would require a resolution of about 10 6 cells for a uniform grid. Clearly, in this case 
either static or adaptive mesh refinement would be useful to resolve the spike. With enough resolution to resolve the 
spike, our algorithms provide an accurate solution in this regime. 

6.3. Radiative shocks with a variable Eddington tensor 

Calculating the structure of radiation modified shocks without i nvoking the Eddington approximation or other 
simplifying assumptions is extremely challenging. ISincell et al.1 ()1999l ) have used a time dependent radiation hydrody- 
namics code to compute the structure o f radiation m o dified shocks. We have found that comparison of our solutions 
in this regime to the results reported by ISincell et al.l (|1999l) is a useful code test. 
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X X 

Fig. 12. — Same as figure llOl but for A4 = 5. At this Mach number, a supercritical shock is formed. Inset in the gas temperature panel 
shows details of the Zel'dovich spike. 

In many ways, the numerical algorithm used by ISincell et al.l (jl999f ) is similar to that described here. They solve 
the radiation moment equations using a variable Eddington factor computed from a formal solution of the transfer 
equation, assuming LTE and gr ay opacities. Howev er, the underlying hydrodynamic solver they use is quite different 
from that adopted in this work. Sinccll et al. (1999) use methods based on artificial viscosity for shock capturing, and 
solve the internal rather than total gas energy equation. The Godunov methods used here use Riemann solvers rather 
than artificial viscosity for shock capturing, and are based on the total gas energy equation. 

Since the Mach number A4 is not known independently of the shock structure, it is not possible to compute the 
solution in the frame of the shock. Instead, we initialize a flow that generates a shock, and allow the shock to propagate 
a large enough distance to settle into a steady structure. To generate the shock, we collide two symmetric flows at 
the center of the domain, producing two identical shocks that propagate symmetrically away from the center. This 
avoids having to implem ent reflecting boundary conditions in the radiative transfer solver that is used to compute the 
VET dDavis et al.ll2012h necessary if the shock is generated by reflection of the flow off a wall. The initial conditions 
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Fig. 13. — Same as figure ITUl but for JVl = 50. At this Mach number, the downstream flow is radiation pressure dominated. 

for the flow are chosen to match those used by iSincell et a l. (1999) (see also SS10) for the subcritical shock solution 
discussed in §2.1 of their paper. We use a computational domain that spans —2 < x < 2. Initially, the density is one, 
the gas pressure is P — 0.923pT, and the gas temperature is T = 1 + 7.b\x\/2. For x < 0, the flow velocity is v = 20, 
while for x > it is v = —20. The dimensionless pressure and speed of light arc P = 1.08 x 10~ 10 and C = 10 6 , and 
7 = 5/3. The radiation temperature is the same as gas temperature, the radiation flux is F r = —(dE r /dx)/(3a a ) 
(the Eddington approximation is assumed in the initial condition only because the density is uniform initially) , where 
the absorption opacity a a = 10. The VET is initialized with only non-zero diagonal components of 1/3; thereafter the 
VET is computed self-consistently with the flow using short characteristics. 

Profiles of various quantities at t = 0.03 are shown in Figure [15] in the region < x < 0.5. (We have checked, as 
another test of our code, that the solution is exactly symmetric with respect to x = 0.) Profiles of each variable show 
the characteristic structure of a subcritical shock, including a strong radiative precursor and Zel'dovich spike. Perhaps 
of greatest interest is the profile of the x — x component of the VET shown in the lower right panel. Upstream of the 
shock f xx is much bigger than 1/3, as expected. Downstream of the shock, it is close to 1/3. However, in the vicinity 
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Fig. 14. — Gas temperature profiles with different Mach number for radiation modified shocks. Blue lines are the initial conditions, while 
black dots are the numerical solution after a few flow crossing times. Note the subtle changes in the shock structure between M = 1 and 
2, which are clearly captured in the numerical solutions. 



of the shock itself, f xx is smaller than 1/3, in agr eement with the res ults of ISincell et al.1 (|1999f ) (see their figure 8). 
The reason for the behavior is discussed in §2.3 of ISincell et al.l (|1999l l. Most of the upstream radiation is generated 
within one optical depth of the shock front. Rays parallel to the shock front go through a longer path of source than 
rays perpendicular to the shock front. Thus, the intensity along rays parallel to the shock front is larger than those 
perp endicular, leading to a value of f xx less than 1/3. Note that since the empirical relations for most flux limiters 
(e.g.. lLevermore fc PomraninsJll98~ll ) bound the Eddington factor between 1/3 and one along the direction of radiation 
energy density gradient, FLD cannot get the right answer in this case. 

We find the most quantitative comparison to the previously published results of lSincell et al.l (|1999h is to plot the gas 
temperature, pressure, and radiative flux against the inverse compression ratio y = Pn/p, where po is upstream density. 
Figure IT51 shows our result at t = 0.3, for d irect comparison to figures 2 and 3 in ISincell et al.l (|1999[). An appr o ximat e 
analytic expression for these quantities bv lZel'Dovich fc Raizerl (I1967I) (and also given in §1.3 of lSincell et all (Il999h l 
is shown as solid lines. The numerical solution of lSincell et al.1 (|1999l ) is in fact in poor agreement with the analytic 
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Fig. 15. — Structure of a subcritical shock computed with a VET at time t = 0.03 in one dimension with 1024 cells. In the upper 
left panel, the red line is the radiation temperature T r while the black line is the gas temperature T. Note that x — x component of the 
Eddington tensor (lower right panel) can be smaller than 1/3 near the shock front. 



predictions, and these authors argue this is due to approximations made in the latter. It is clear from Figure [T6l that 
our solutions agree very well with the approximate analytic solution, and therefore do not agree with the results of 
iSincell et al.l ()1999t ) for these quantities. Recall these authors use a numerical method based on artificial viscosity for 
shock capturing, and the internal rather than total energy equation. The plot measures the internal structure of the 
shock, which could be quite sensitive to both the form of viscosity used, and the degree to which energy is conserved. 
Our solution may be in better agreement with the analytic expectations since our algorithm uses a Riemann solver 
rather than artificial viscosity to capture shocks, and is based on the total energy equation. 

6.4. Photon Bubble Instability 

Photon bubble instability is an overstable mode present in radiation sup ported, magnetized a tmospheres. It was first 
noticed in numerical simulations of neutron star polar cap accretion flo ws (iKlein fc Aronsl[l989i) , and was subsequently 
confirmed in a linear stability analysis for neutron star atmospheres (|Aronslll992[ ). More recently, iGammiel ([19981 ) 
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Fig. 16. — Gas temperature, radiation flux, and gas pressu re a s a function of inverse compression ratio r\ = po/p, where po is upstream 
density, for the subcritical radiative shock shown in Figure 1151 at time 0.3. Temperature and pressure are scaled by the downstream 
temperature T\ and pressure Pi at x = 0, while the radi ative flux is scaled b y the flux F r , a at the shock front. Note F r is the co-moving 
flux. The red lines are analytic solutions given bv lZel'Dovich fc Raiz cr (1967) (see also lSincell et al.|[l999h . 



and iBlaes fc Socrates! (|2003h recognized a similar instability occurs in accretion disk atmospheres, and analyzed the 
properties in the li near regime. Numerical simulations that in vestigate the non-linea r regime were performed by 
iTurner et all (|2005| ) using the FLD module in the ZEUS code (|Turner fc Stone! I2001D . Since both the linear and 
non-linear regimes of this instability have been studied in detail in the literature, and since fundamentally it relies on 
radiation modified MHD waves, it provides is an excellent test of our radiatio n MHD algorithms in multidimensions. 
In this section, we repeat one set of parameters reported bv lTurner et al.l (|2005l ) and compare our results with solutions 
from ZEUS. In order to make direct the comparisons, we adopt the Eddington approximation, so that the diagonal 
components of the VET are fixed to be 1/3. 

The initial conditions consist of a stratified atmosphere in mechanical and thermal equilibrium, in which a uniform 
radiation flux balances gravity. The structure of the atmosphere is computed from the solution of the hydrostatic 
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equilibrium equations 

^7- = -pg + vFr, P = pT, 
dy 

^ = -*F r , E r — T 4 . (27) 
3 dy 

with a constant gravitational acceleration g = 1 in the — y direction, and an opacity of a = 137.11/3+ 0.15p 2 T~ 3 5 
(which includes both electron scattering opacity and Kramer's opacity). The co-moving vertical radiation flux F r ,o = 
2.52 x 1CP 4 . These equations are integrated in a domain of size (0,2.75) x (26.11,28.86) in our dimensionless units, 
starting from the midplane y = 27.485 where the dimensionless pressure, temperature, density and radiation energy 
density are all set to 1. In this test, the dimensionless parameter P = 762 (so the ratio between radiation pressure and 
gas pressure in the midplane is 254), and the dimensionless speed of light is C — 7557.74. The initial conditions include 
a uniform, horizontal magnetic field in the +x direction with magnetic pressure 10% of mid-plane radiation pressure. 
Our simulations use a resolution of 256 x 256 grid points. The instability is seeded with random perturbations of the 
density with amplitude 10~ 3 at the mid plane. The initial condition for this simulation is similar (but not identical, 
the am plitude of perturbation is 10~ 8 in lTurner et al.ll2005f) to the initial condition used for Figure 9 of I Turner et al.l 
(2005). This model is appropriate for the surface layers of an accretion disk at 20 Schwarzschild radii from a 10 8 M Q 
black hole. In this case, the middle of the simulation domain corresponds to a location one scale height above the 
disk midplane, and the simulation domain extends from 0.95H to 1.05-ff. The total optical across the domain in the 
vertical direction is 364, thus the Eddington approximation should be an excellent description. In our dimensionless 
units, one orbital period is 6.28 time units. 

To measure the linear growth rate of the photon bub ble instability as a function of the vertical and horizontal 
wavenumbers k x and k y , we follow the method used by iTurner et all (|2005l ). We Fourier transform the horizontal 
velocity v x at 50 snapshots between times of t — 0.8 and 2.8 equally spaced at a time interval of 0.04. We assume 
exponential growth at each wavenumber independently, and therefore calculate the slope of the logarithm of the power 
with time between each neighboring time, and average the resulting rates from all the snapshots together. This 
approach only approximates the true growth rates in a stratified atmosphere. The result is shown in in Figure I17[ 
which can be com pared directly to Figure 9 of ITurner et al.l ()2005l ) (which also shows solution of linear analysis from 
iBlaes fc Socrates! i|2003h ). There is significant noise in the measurement of the growth rates, especially at small k. 
However, clearly the correct trends predicted by the linear ana l ysis a re reproduced, with the fastest growth occurring 
for modes with k x = k y . The noise in Figure 9 of ITurner et al.l (|2005| ) is much smaller because they can use an initial 
amplitude as small as 10~ 8 . This cannot be achieved in Athena because the pressure gradient and gravity cannot be 
balanced to roundoff error (without special modification to the reconstruction algorithm) and there is always noise in 
the background state. 

In the non-linear regime, ITurner et al.1 (|2005h found the instability produces shock trains that propagate from the 
bottom to the top of the atmosphere. With a purely horizontal background field, there is no preference for shocks 
propagating in either direction, resulting in a symmetric pattern of shocks in the atmosphere. The left panel of Figure 
IT8l shows the density and velocity distribution at time 9.83 when the shock train pattern becomes relatively strong. 
The velocity vectors are nearly horizontal direction because flow is confined to be parallel to the s trong horizontal 
magnetic field. This figure can be compared with the first panel of Figure 19 of ITurner et all (12001 : our results are 
very similar. 

In the FLD approximation used by ITurner et al.1 ([2005), the co-moving radiation flux is always assumed to be 
along the direction of the gradient of the co-moving radiation energy density. Although we adopt the Eddington 
approximation in our simulation, we make no assumption about the direction of the flux. Thus, it is possible for us 
to test whether the radiation flux is indeed everywhere parallel to the gradient of E r in this highly nonlinear flow. 
Figure [18] shows the distribution of cos 6* at t = 9.83, where the angle 8 is the angle between the perturbed co- moving 
radiation flux in our simulation F r ,o — F r fl and the direction of the flux predicted by FLD — V_E r .o/(3cr t ). If the 
FLD approximation is valid, cos8 will be one everywhere. We find that near the shock fronts, where there is a large 
density gradient, this assumption clearly fails. Since th e optical dept h acro ss the domain is so large, the Eddington 
approximation should apply, and the flux-limiter used by ITurner et al.l (|2005f ) should have a numerical value very close 
to the value of 1/3 adopted here. Thus, the discrepancy in the direction of the radiation flux is not due to the value 
of the Eddington factor, but rather because in FLD the radiation flux is always along the gradient of radiation energy 
density. 

6.5. Shadow Test 

Since the direction of the flux is not known independently of the energy density, it is well-known that methods based 
on the diffusion approximation cannot represent shadows. Thi s is often demonstrated through the irradiatio n of very 
optically thick structures with a beamed radiation field (e.g.. iHayes fc Normanl 12001 iGonzalez et al.ll2007f ). In this 
section, we demonstrate our methods capture shadows accurately, and moreover we present the dynamical evolution 
of a cloud ablated by an intense radiation field. VET is crucial to to get t he solution correct ly, which is calculated 
from the transfer module using short characteristic. The tests described in iDavis et al.l (|2012t l (see Figure 6 of that 
paper) show the accuracy of the radiation transfer solver for two radiation beams. Here we show that when coupled 
to the VET method, our solver captures shadows correctly. 




Fig. 18. — Non-linear outcome of Photon Bubble instability at time 9.83, corresponding to 1.56 orbital periods. The left panel shows 
the density and velocity vectors. The right panel plots cos 9, where 9 is the angle between the perturbed co-moving radiation flux and the 
gradient of the perturbed co-moving radiation energy density. Deviations of cos 9 from one near the shock fronts indicate the difference 
between the direction of the flux in our simulations versus the value assumed in FLD. 
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Our test is performed in a rectangular box of size (—0.5,0.5) x (—0.3,0.3) cm. Initially the background medium 
has density po = 1 g/cm 3 and temperature To = 290 K. An over-dense clump is located in an elliptical region 
r = x 2 /a 2 + y 2 /b 2 < 1, with a — 0.1 cm and b = 0.06 cm. The density inside this region is p(x,y) = po + (pi — 
Po)/[1-0 + exp(10(r — 1))] with p\ — 10po- The clump is in pressure equilibrium with its surroundings, so the interior 
is much colder than the ambient medium. The initial radiation temperature is the same as the gas temperature 
everywhere. Pure absorption opacity is used with a — (T/To) _3 ' 5 (p/po) 2 cm -1 . The radiation flux F r is zero 
everywhere initially. Outflow boundary conditions are used on both y boundaries, and on the right x boundary. A 
constant radiation field with temperature T r = 6Tq is input through the left x boundary at angles of ±14 degrees with 
respect to the x axis. This radiation field is also input along the top y boundary at —14 degrees, and along the bottom 
y boundary at +14 degrees, in order to mimic the radiation field from an infinite plane along the left x boundary at 
x = —0.5. At the left x boundary, the gas temperature and density are fixed to To and po respectively. A resolution of 
512 x 256 cells is used for this test. The dimensionless speed of light C = 1.9 x 10 5 , and the parameter P = 2.2 x 10 -15 . 

In the ambient medium, the photon mean free path is the length of simulation domain, so that the radiation can 
propagate freely across the box. However, inside the clump the photon mean free path is only 3.2 x 10~ 6 of the length 
of simulation box due to the high density and low temperature. Thus, radiation coming from the left boundary cannot 
penetrate the clump, and shadows are cast towards the right boundary. By using an input radiation field at two angles, 
both umbra and penumbra are formed. This makes the test more difficult, because ad-hoc closures that capture only 
one direction for the flux will not represent both the umbra and penumbra correctly. 

In Figure [19] , we show the radiation energy density, and x — x and y — y components of the VET after one timestep 
at St = 1.0e~ 3 . This timestep corresponds to the CFL stability condition set by the adiabatic sound speed in the 
ambient gas and the grid resolution. The ratio of this timestep to the light crossing time of the domain is about 500, 
thus after one timestep both radiation beams have crossed the domain and a steady-state shadow structure is formed. 
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Fig. 19. — Shadows created by the irradiation of an optically thick clump by two beams at ±14 degrees with respect to the x axis. From 
top to bottom, the three panels show the radiation energy density E r , and the xx and yy components of Eddington tensor f xx and f yy 
respectively. Inside the clump, the Eddington tensor is close to 1/3 due to the large opacity, while E r retains its initial value. Umbra and 
penumbra are clearly formed behind the clump. 



Since the radiation temperature is low in this first test, the clump undergoes negligible dynamical evolution. To 
demonstrate the dynamics of the ablation of a clump, we have repeated this test with an increase the background 
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temperature To by a factor of 357 such that the dimensionless P = 10~ 7 . Then the radiation pressure of the applied 
radiation field at the left x boundary is about 10~ 3 of the gas pressure. This non-negligible radiation field pushes and 
compresses the clump from the left and the shape of the clump is changed. A shock is driven into the clump during 
this process. The clump is also heated up by the radiation field and outflow from the surface of the clump is formed on 
the left side. The density distribution and velocity field at three different times are shown in Figure [20] The adiabatic 
sound speed in the ambient medium is 1.3 in our units and 0.39 inside the clump. In the first panel of Figure [20} the 
density on the left side of the clump is increased by a factor of 2 due to the radiation pressure. Outflow is generated 
from the surface of the clump as shown by the velocity vectors. The flow is subsonic in the ambient medium and 
supersonic with respect to the sound speed inside the clump. In the second panel, we clearly see a shock driven to 
the clump with Mach number up to 3. In the third panel, the clump is completely destroyed by the radiation field 
and several shock fronts are forme d. The ablation of clouds by strong radiation field is in fact of great interest in 
the interstellar medium (ISM) (e.g.. iBertoldi fc McKedfl99(l lBallvlll995| ) and near the central region of quasars (e.g., 
lMathewsl lT983) . Further explore of this flow is warranted, but beyond the scope of this paper. 
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Fig. 20. — Ablation of an optically thick clump by two radiation beams at ±14 degrees with respect to the x axis. From top to bottom, 
the density (color) and velocity (vectors) are shown times t = 4.9 X 10~ 7 , 1.3 X 10~° and 2.6 X 10 — 6 s respectively. 



7. DISCUSSION 

7.1. Comparison with methods based on the diffusion approximation 

A central ingredient of the algorithm presented in this paper is that it uses a VET to close the radiation moment 
equations, rather than relying on the diffusion approximation. It is therefore worthwhile to discuss specific examples 
where there will be differences between these approaches. 

The first difference is that by solving the moment equations for the radiation flux, we are able to solve for the 
direction of the flux self-consistently with the flow, rather than assuming it lies in the direction of Vi? r . Both the 
photon bubble instability (see Figure fTHT) and the shadow test (see Figure [T9]) demonstrate the importance of keeping 
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the direction of the flux. 

A second difference is that the value of the flux limiter is always an approximation (of unknown reliability) to the 
true value of the diagonal components of the VET. In fact, the planar subcritical shock test (see Figure [T5t shows 
that near the shock front, the VET can be less than 1/3 along the direction of the r adiation energy density gradient , 
which clearly cannot be represented by the formulae used for most limiters (e.g., iLevermore fe Pom raning 19811). 
Based on the density dist ribution from a snapshot of simulations for black hole accretion disk done by iHirose et al.1 
(200|), iDavis et al.1 (12011) compares the Eddington tensor used in the flux-limited diffusion approximation and the 
VET calculated based on our transfer module. They also show significant difference (see Figure 7 of that paper). In 
the photosphere, the Eddington tensor given by the diffusion approximation can be much larger than the values given 
by the VET. 

A third difference is that we do not drop the time-derivative of the flux dF r /dt in the moment equations for the flux. 
This term is important when the inertia of the radiation field is important, and dropping this term can lead to incorrect 
behavior in this regime. For example, consider a uniform medium with density p = 1, absorption opacity a a = 10 
and scattering opacity a s — 10 in thermal equilibrium (E r = X" 4 = 1 in our dimensionless units) with zero fixed fram 
radiation flux, and an initially uniform velocity along the x-direction vq = 1. Adopt values for the dimensionless speed 
of light C = 100 and parameter P = 1000, and use the Eddington approximation. In the diffusion approximation, 
because the radiation energy density is always uniform, the co-moving radiation flux will be always be zero and the 
system will always be in thermal equilibrium. However, this is not true if the dF r /dt term is kept. Since the gas has 
a non-zero velocity with respect to the radiation field, the co-moving radiation flux is not zero, and there will be a 
drag force from the radiation field on the fluid, (a consequence of the fact that solutions to the equations of radiation 
hydrodynamic are not Galilean invariant). Because of the drag with the radiation field, momentum will be exchanged 
and the gas will be heated. In this example, this process can be described by 

9 (H_ ro / F? vE r + v-P T 



dt C 

^_ Cw ( Fr _*±^). (28) 

In solutions to these equations, p and E r are unchanged while the change of velocity with time t can be described by 
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We compare this analytic solution for this problem, and a numerical solution generated by our algorithms, in Figure [2"T1 
Our solutions captures this momentum exchange process accurately, which is not possible with flux-limited diffusion. 
The importance of this process is characterized by the parameter P/C 2 , which measures the importance of momentum 
carried by the radiation field compared with the material momentum. 

A fourth difference between our VET algorithm and FLD is that we keep the non-diagonal components of the 
VET. In the FLD approximation, th e radiation pressure is always a diagonal tensor. Thus, radiation viscosity (e.g., 
lMihalas"fc Mihalas 1984; Castor 2004) cannot be captured in this approximation. However, with the VET approach, 
the off-diagonal components are computed self-consistently with the flow. Thus, the effects of radiation viscosity are 
captured naturally in simulations of accretion flows. 

The importance of the differences between VET and FLD depend on the particular application of interest. For 
example, in the inner regions of accretion disks around supermassive black holes, we expect the momentum carried by 
the radiation field t o be a significant frac tion of material momentum (P/C 2 is not negligible), and photon viscosity may 
be significant (e.g,. IAgol fc K rolik 19981). Thus, adopting algorithms based on VET is important for our applications. 

7.2. Performance and Scaling 

For solving the application problems of interest, our algorithms must not only be accurate, but also efficient. In this 
subsection, we report preliminary performance and scaling of the method on parallel systems. We continue to work 
on optimization of the algorithms, and it is likely substantial improvements in performance are still possible. 

For reference, it is useful to compare to the performance of Athena in ideal MHD without radiation. On a single core 
of a 2.5 GHz Intel Xcon processor, Athena updates roughly 1 x 10 5 grid cells per second for 3D MHD calculations with 
a 32 3 grid. Weak scaling is very good, with an efficiency (defined as the performance per core in a parallel calculation, 
compared to the performance of a serial calculation) of 0.9 on up to 10 5 processors. 

With radiation hydrodynamics and MHD, there are two potential bottlenecks. The first is the implicit solution of the 
linear radiation subsystem, which requires inverting one matrix per timestep (this is already an advantage over methods 
that require implicit solution of nonlinear systems which generally require many matrix inv ersions per tim estep) . To 
inver t the matrix, we hav e explored the use of algorithms implemented in both the LIS (e.g ., [N ishida 2 0l"ol) a nd Hypre 
(e.g.- lFalgout et aDl2006f l libraries, as well as our own implementation of the multigrid (e.g.. lPress et al.ll"l992[ ) method. 
Both libraries provide a variety of iterative algorithms and preconditioners for solving sparse matrices, and can be 
used on both serial and parallel systems. The convergence rate of the matrix solver is limited by several things, such 
as the opacity and error tolerance. 

To check the serial performance of each of these methods, we use the linear hydrodynamic wave problem as described 
in Section UTTl on a 3D domain and a 32 3 grid. With the LIS library on a single 2.5 GHz Intel Xeon processor with 
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Fig. 21. — Momentum exchange between radiation and material due to radiation drag. Red line is the analytic solution given by equation 
1291 while the black lines are our numerical solutions for normalized velocity and radiation flux. The solution reaches equilibrium when the 
co-moving radiation flux is zero. The velocity changes by 12% in this example. 



tolerance level 10 -8 , we achieve roughly 1.7 x 10 4 cell updates per cpu second in optical thick regimes with the ILU 
pre-conditioner. The Hypre library is somewhat faster, achieving about 3 x 10 4 cell updates per cpu with the GMRES 
solver and the BoomerAMG preconditioner. Because the sparsity pattern of the matrix we invert is not common (see 
Appendix A), we must use the general matrix solvers in these libraries, which have very poor cache performance. By 
writing our own multigrid algorithm for inverting the matrix, we can take advantage of the specific structure of our 
matrix to substantially improve performance. With multigrid, we can achieve 4.3 x 10 cell updates per cpu for this 
test, and the convergence rate is independent of whether the linear waves are optically thick or thin. With both the 
LIS and Hypre libraries, performance is about 2^5 times slower in the optically thin regime. Currently, the multigrid 
solver is our workhorse algorithm. 

A second bottleneck to performance in simulations with radiation is the formal solution of the transfer equation. 
The cost of this step depends on the number of rays required to resolve the angular dependence of the radiation field, 
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and whether the LTE approximation can be used to compute the source function fsee lDavis et al.ir2012f) . In our tests, 
we have found sufficient accuracy is obtained with roughly 24 angles in 3D, although this number may be highly 
application dependent. If absorption opacity is dominant and the LTE approximation can be used, this step consumes 
a negligible fraction of the overall execution time. Moreover, since the cost to compute the VET scales linearly with the 
number of angles, using many more angles (for example, 80 instead of 24) does not affect the overall cost per timestep 
significantly. For example, using 80 angles increases the cost per timestep by less than factor of two compared to 24 
angles. 

However, if scattering opacity is dominant and the non-LTE integrator must be used to compute the VET, then 
this step can dominate the cost. With 24 angles in a non-LTE calculation, the overall performance of the code can be 
decreased by a factor of up to 5. Moreover, using more angles decreases the convergence rate of the iteration scheme 
used for non-LTE problems, thus roughly tripling the number of angles (from 24 to 80) will increase the cost per 
timestep by more than a factor of three. There are several simple steps that can be used to increase performance in 
this regime, however. For example, in many cases it is possible to compute the VET only every several MHD time 
steps, since the geometry of the radiation field (and therefore the VET) changes only slowly compared to the flow 
structure. 

It is interesting to compare the performance of the VET algorithm to FLD. Since we have not implemented FLD 
directly in Athena, we cannot perform a head-to-head comparison. However, w e note that because the implicit 
differencing most often used with FLD (e.g. jTurner fc S tone 200ll; IZhang et ai1l2011l ) involves nonlinear terms, Newton- 
Raphson iterations are required with one matrix solver per iteration. Depen ding on the convergen ce rate of the 
iterations, this could require several to many matrix solvers per timestep. Thus. [Turner fc Stond (|2001h reported that 
in the optically thin limit the FLD module in ZEUS made the code run 10 x slower, whereas we find the VET algorithm 
(including the short-characteristics RT solver) is only 2 — 3x slower than the standard MHD integrators in Athena. 
ZEUS itself is several times faster than Athena, so the VET may not in fact be faster than FLD, but it is at least 
comparable in performance. Using FLD because it is faster does not seem to be supported by these comparisons. 

The parallel performance of the code is also limited by the matrix solver we use. Generally, in parallel calculations, 
the iterative linear system solver requires more steps to converge compared to the serial case. Using the linear 
hydrodynamic wave problem as described above, with 32 3 grids on each processor, we have performed a weak scaling 
test for both the LIS and Hypre libraries, and our own multigrid solver. With 8 cores, the efficiency of the two libraries 
are about the same: 0.65, and on 256 cores this drops to 0.25. With our multigrid solver, the efficiency on 8 cores 
is 0.76, and drops to 0.56 on 256 cores. While this is substantially lower than the efficiency of the MHD integrator 
alone, it is still sufficient to enable large 3D applications. However, clearly further effort to improve performance is 
warranted. 

8. SUMMARY 

Wc have described a Godunov algorithm for multidimensional radiation hydrodynamics and MHD, designed to study 
the dynamics of radiation dominated accretion flows around compact objects, and rad iation driven wind s from stars, 
disks, and galaxies. The algorithm has been implemented in the Athena MHD code (|Stone et al.ll2008h . and tested 
over a wide range of parameter space. 

Compared with previous algorithms for radiation hydrodynamics, there are several features of our approach that 
bear special mention. 

Use of a VET rather than FLD. Most multidimensional algorithms for radiation hydrodynamics adopt the FLD 
approximation. However, this means that the direction of the radiation flux is lost, that the value of the flux is 
determined by an ad-hoc limiter, that the inertia of the radiation field is lost and that the off-diagonal components of 
the radiation pressure are assumed to be zero. We have shown through a variety of test problems that using methods 
based on a VET overcomes these limitations. 

Formal solution of the transfer equation via short characteristics. To compute the VET, we discretize the specific 
intensity over many angles, and solve the transfer equation using short characteristics. While this has been deemed 
prohibitively expensive for 3D applications, advances in algorithms and hardware have overcome this difficulty. We 
find we can solve the transfer equation in full 3D using 24 angles at a cost that is no more than a single timestep of 
the MHD algorithm. We have also explored the use of Monte Carlo methods to compute the VET, and have found 
for our applications, short characteristics is substantially faster (| Davis "etld1l27)l2h . 

Use of Godun ov methods for th e underlying MHD solve r. The VET m ethod has been implement ed previously in 
the ZEUS (e.g.. IStone et alj|1992t fHaves fc Normanll2003l ). and TITAN (jGehmevr fc Mihalasl [l99l codes, however 
these methods used an operator split method with an artificial viscosity for shock capturing. Instead, we have adopted 
a higher-order Godunov method with a dimensionally unsplit integrator. Comparison of the structure of radiation 
modified shocks show the advantages of this approach. Since the addition of stiff source terms to Godunov methods 
is problematic, we extend the modified Godunov method of SS10 to multi-dimensions. We show this method results 
in accurate and stable evolution. 

To test our algorithms, we have introduced a variety of quantitative test problems. We show the results of error 
convergence tests for linear radiation modified acoustic waves in both hydrodynamics and MHD over a wide range 
of parameters, and give examples of eigenmodes that can be used by others for testing. We also use the structur e 
of radiating shocks in non-equilibrium diffusion computed using the methods outlined bv lLowrie fc Edwar ds (2008). 
Our algorithm reproduces the subtle changes in shock structure that occur as the Mach number is increased, even 
including very high Mach numbers where radiation pressure dominates the postshock flow. We also have performed 
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multidimensional tests based on the photon bubble instability and shadowing. 

There are still some limitations to our algorithms, which could be improved with future development. For example, we 
assume LTE, and our methods are limited to the grey approximation. Extension to multiple frequency bands requires 
analytic expressions for the interaction between different frequency bands in order to calculate the propagation operator 
used to update the stiff source terms. This is an important direction for future work. Our algorithm cannot propagate 
linear waves very well when absorption opacity is very large, in particular when r a > C/P. Note that the algorithm is 
still capable of capturing flows in this regime. However, the damping rate and phase velocities of radiation modified 
acoustic waves will not be accurate. This can be improved by solving the total conservative form of the equations 
in this regime. If absorption opacity is so large that radiation temperature is locked to the gas temperature, the 
cr a (T 4 — E r ) term can be dropped in the radiation energy source term. Then damping rate and phase velocity of the 
radiation modified acoustic waves can be captured accurately in this regime. Since we solve the radiation transfer 
equations in the mixed- frame, our algorithm cannot be used for problems when spectral lines dominate the transport. 
The current algorithm is also not appropriate for problems when the radiation field changes rapidly (on a light-crossing 
time) because the VET is calculated using the time-independent transfer equation, although the error introduced by 
this assumption is acceptable if the time step is comparable to the light-crossing time of each cell, as we demonstrated 
in Section [OJ 

Not all the applications of interest requiring including the effects of radiation pressure. In this regime, it is possible 
to use a much simpler algorithm which only adds the radiation source terms to the energy eq uations, computed from 
a formal solution of the transfer equation. This method is described in a companion paper (|Davis et alJl2012t l. and 
is currently being used to study black hole accretion disks in the gas pressure dominated regime. We use the same 
module for computing the VET in this work. 

Currently we are using the radiation MHD version of Athena to study several applications, including the nonlinear 
regime of the Ray leigh- Taylor instability in a radiation supported atmosphere, radiation feedback in AGN, and the 
saturation of the magneto-rotational instability (MM) in radiation dominated disks (Jiang et al. in preparation). 
Eventually the code will be publicly available through the Athena Trac site. 
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Implicit backward Euler differencing of the radiation subsystem (equation [20] and [2Tj) results in a set of linear 
equations for the radiation energy density and flux in each cell. For every cell with indices (i,j,k), we have the 
following four equations 
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APPENDIX 



MATRIX FOR THE 3D RADIATION SUBSYSTEM 



e E? +1 (i,j,k- l) + 9 1 F?+ 1 (i,j,k- 1) 

+ e 2 E?+\i,j - i,k) + e 3 F? 2 +1 (i,j - i,fc) 

+ 9 4 E? +1 (i - 1 , j, k) + 9 5 F^ +1 (i-l,j,k) 
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cj> E? +1 (i,3, k - 1) + cPiFf+^iJ, k - 1) 
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+ *Pi2E? +1 (i, j, k + l) + ^F? 2 +1 (i, j, k + 1) 



tp E? +1 (i,j, k-l) + Vl F^ l (i,j,k- 1) 



(A4) 



+ ^E^{i,j -l,k) + V 3 F^ l (i,j -l,k) 
+ ^ +1 (i - l,j,k) + ^F?+\i -1,3, k) 
+ WE?+ I (i,3,k) + IP7 F?+ I (i,3,k) 
+ p 8 E>; +1 (i + l, j, k) + ^F^ 1 (i + l,j, k) 
+ Vw El l+1 {i,j + l,k) + <pnF? 3 +1 (i,j + l,k) 
+ Vi2K +1 {i, 3, k + l) + <p 13 F? 3 +1 (i, 3, k + 1) 
= F? 3 (i,3,k)+d 2 a a v z f^ k /C. 

Here v x , v y and v z are the velocities along the x—, y— and z~ directions respectively, and d 2 = AtC. The temperature 
T[j k is the value we estimate in Section 13.3.11 The velocities and opacity are at time step n and are kept constant 
during the solution of the radiation moment equations. 

Let the components of the symmetric 3x3 VET f be labelled as 



f : 



/ll /l2 /l3 
/21 fl2 /23 
/31 /32 /33 



We also define the following quantities 

d lx = AtC/{2Ax), d ly = AtC/{2Ay), d u 



CM: 



C301 



CkO: 



ftt 2 (i,j,k)-fU 2 (i-l,j,fc) 
fli\i,3,k) + fl[\i-l,j,k) 
llj 2 {i,3, k) - f 22 2 (i,j - 1, k) 
fll 2 {h3\ k) + fU 2 (i,j - 1, k) 
fU 2 ^,3,k)-fU 2 {z,3,k-l) 
ill 2 {hi, k) + fH 2 {i,j, k-1) 



Cil 



Cjl 



Ckl 



= A«C/(2Az), 
fU 2 (i + l,j,k)-fU 2 (i 



(A5) 



(A6) 



fll 2 {i+l,3,k) + fl[\i 



fll\i,3,k + l)-fli 2 {i,j,k) 



3,k) 



3,k) 
,3,k) 



J33 (*>J' fc + 1 ) + /33 ih3\k) 
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Then the coefficients 9q — #15, (f>o — tfii 3 , tfio — cf>i 3 and cpo — ipi 3 in equations IA1I — PQI are 

9 Q = -d lz (l + CkO)f 3 ( 2 (i,j,k-l), 

91 = -d lz (l + CkO), 

9 2 = -d ly (l + CjO)f^ 2 (i lJ -l,k), 

9 3 = -d ly (l + CjO), 

0i = -d lx (l + CiO)fU 2 (i - l,j,k), 
6 B = -d lx (l + CiO), 

e 6 = i + d lx (i + ai)fl( 2 (i,j, k) + d lx (i - ao)fl( 2 (i,j, k) + d ly (i + Cji)fV%j, k) 

+ d ly (l - CjO)fU\i,j, k) + d u {l + Ckl)fU%j, k) + d u (l ~ CkO)fU\i,j, k) 
+ Kd 2 ((T a - a s )v x [v x (l + fu(i,j, k)) + v y f 12 (i,j,k) + v z f 13 (i,j,k)} /C 2 
+ Kd 2 ((T a - a s )v y [v v (l + f 22 (i,j,k)) + v x f 12 (i,j,k) + v z f 23 (i,j,k)} /C 2 
+ Kd 2 {a a - cr s )v z [v z (l + f 33 (i,j,k)) + v x f 13 (i,j,k) + v y f 23 (i,j,k)} /C 2 
+ lZd 2 a a , 

7 = d lx (l + Cil) - d lx (l - CiO) - 1Zd 2 v x (o a - a.)/C, 
9 s = d ly (l + Cjl) - d ly (l - CjO) - Kd 2 v y (a a - a B )/C, 
9 9 = d lz (l + Ckl) - di,(l - CkO) - Tld 2 v z {a a - a.)/C, 

e w = -d lx (i-ca)fl( 2 (i + i,j,k), 

011 = ^(1 -Cil), 

e 12 = -d ly (i-cji)fU 2 (i,j + i,k), 

6 13 = d ly (l-Cjl), 

On = -d u (l - Ckl)fU 2 (iJ, k + 1), 
9 15 = d lz {l-Ckl); 

cf> = -d lz (l + CkO)f 31 (i,j,k-l), 

cf>i = -d lz (l + CkO)fU 2 (i,3,k-l), 
<h = -di v (l + CjO)f 2 i(i,j-l,k), 

c(> 3 = -d ly (l + CjO)f^ 2 (i,j -l,k), 
(t> 4 = -d lx (l + CiO)f 11 (i-l,j,k), 

<k = -di*(l + CiO)fU 2 (i-l,j,k), 

fa = d la (l + Cil)fi\{i,j, k) - d lx (l - CiO)f n (i,j, k) + d ly (l + Cjl)f 21 (i,j, k) 
-d ly (l - Cj0)f 21 (i,j,k) + d lz (l + Ckl)f 31 (i,j,k)-d lz (l - Ck0)f 31 (i,j,k) 
-d 2 a t [v x (l + fu(i,j, k)) + v y f 2 i{i,j, k) + v z f 3 i(i,j, k)] /C + d 2 a a v x /C, 

<h = 1 + di x {l + Cil)fll 2 (i,j, k) + d lx (l - CiO)fU 2 (i,j, k) 

+ d ly (l + Cjl)fU 2 (i,j, k) + d ly (l - CjO)fU 2 (i,j, k) + d u (l + Ckl)fU 2 {i,j, k) 

+ d u (l - Ck0)f 3 ( 2 {i,j, k) + d 2 a t , 
tfi s = di x (l - Cil)fn(i + l,j,k), 

cj> 9 = -d lx (l-Cil)fll 2 (i + l,j,k), 
4>io = di y (l - Cjl)f 21 (i,j + l,k), 

<Pu = -d ly (l - Cjl)fy 2 (i,j + l,k), 
</>i 2 = d lz (l-Ckl)f 31 (i,j,k + l), 

<ha = -duO- ~ Ckl)fU 2 (i,j,k + 1); 
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^ = -d lz {l + Ck0)f 32 (i,j,k- 1), (A9) 

Vi = + CMi)fU\i,3, k-1), 

4> 2 = -d ly (l + CjO)f 22 (i,j -l,k), 

1 p 3 = -d ly (l + CjO)f 2 1 ( 2 (i,j -l,k), 
i>4 = -d lx (l + CiO)f 21 (i - k), 

i> 5 = -d lx {l + CiO)fl( 2 (i-l,j,k), 

ip 6 = d lx (l + Cil)f 21 {i,j,k) - d lx (l - CiO)f 21 {i,j,k) + d ly (l + Cjl)f 22 (i,j,k) 
- d ly (l - Cj0)f 22 (i,j, k) + d lz {\ + Ckl)f 32 (i,j, k) - d lz (l - Ck0)f 32 (i,j, k) 
-d 2 a t [v y (l + f 22 (i,j, k)) + v x fi 2 (i,j, k) + v z f 32 (i,j, k)] /C + d 2 a a v y /C, 

i> 7 = 1 + d lx (l + Cil)fl( 2 (i,j, k) + d lx (l - CiO)fl( 2 (i,j, k) 

+ d lv (l + Cjl)fU 2 (i,j,k) + d ly (l - CjO)fU 2 (i,j,k) 

+ d u (l + Ckl)fli 2 {i,j, k) + d lz (l - CkO)fU 2 (i,j, k) + d 2 a u 
i> s = d lx (l-Cil)f 12 (i + l,j,k), 

i> 9 = -d lx (l-Cil)fl( 2 (i + l,j,k), 
^io = rfij,(l - Cjl)f 22 (i,j + l,k), 

<M = -di„(i - Cji)fH 2 {i, j + l, k), 
il>i2 = di z (l — Ckl)f 32 (i,j, k + 1), 

Via = - Ckl)fU*(fi,j, k + l); 



Vo = -d u (l + Ck0)f 33 (i,j, k-1), (A10) 

¥'2 = -rfi J/ (l + CiO)/ 23 (i,j-l,fe), 
¥33 = -rfi ?/ (l + CiO)/ 2 1 2 /2 (i,i-l,fc) J 
954 = -^ix(l + Ci0)/ 3 i(i - k), 
ip5 = -di x (l + CiO)ft( 2 (i-l,j,k), 

ipe = d lx (l + Cil)f 31 (i,j, k) - d lx (l - Ci0)f 31 (i,j, k) + d ly (l + Cjl)f 32 (i, j, k) 
-d ly (l - CjO)f 32 {i,j,k) + d lz (l + Ckl)f 33 (i,j,k) - d lz (l - Ck0)f 33 (i,j,k) 
-d 2 cr t Ml + /33(i, j, k)) + v x f 31 (i,j, k) + v y f 32 (i,j, k)] /C + d 2 a a v z /C, 

<p T = l + d lx (l + Cil)fl( 2 (i,j, k) + d lx {\ - CiQ)fl( 2 (i,j, k) 

+ d ly (l + Cjl)fU 2 (i,j,k) + d ly (l - CjO)fU 2 (i,j,k) 

+ di»(l + Ckl)fU%j, k) + d lz (l - CkO)fU 2 (i,3, k) + d 2 a u 
¥>8 = di x (l - Cil)f 31 (i + l,j,k), 

ip 9 = -d lx (l-Cil)ft( 2 (i + l,j,k), 
<Pio = d ly (l - Cjl)f 32 (i,j + l,k), 

<pii = -di y (l-Cjl)f 2 £ (i,j + l,k), 
<Pi2 = d lz (l - Ckl)f 33 (i,j,k+ 1), 

Vi3 = -dix(l ~ Ckl)fU 2 (i,j, k + l). 
The above equations can be written in matrix form, with the vector of unknowns containing elements 
(E r , F r \, F r2 , F r3 )"~j k with index i varying fastest. The result is a large sparse banded matrix that must be in- 
verted (see section I7T2)) . Note that some of the elements of the matrix depend on the boundary conditions. We have 
implemented periodic, reflecting, inflow and outflow conditions. 

DISPERSION RELATION FOR LINEAR WAVES IN RADIATION HYDRODYNAMICS 

In this appendix, we derive the dispersion relation for linear waves in radiation hydrodynamics used for the tests 
in section T6. II The opacity a a — &t = & is a constant. The background state is chosen to be po = l,Po = L-Ev,o — 
l,Xo = 1, which is uniform with zero velocity and radiation flux, and is in mechanical and thermal equilibrium. 
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Following the analysis done by I Johnson fc Kleinl (|2010t ). for each perturbed quantity, we look for wave solutions of 
the form Sf = Sf x e ^t-kx)^ where fc is the wave number, oj the angular frequency and Sf the amplitude. Note that 
we only need the real part of Sf. The amplitude Sf can be a complex number, the imaginary part of which represents 
the relative phase shift between different quantities. Then the dispersion relation between oj and fc has the form 

c 4 fc 4 + c 2 k 2 + c = 0, (Bl) 

where the coefficients are 

7 ojC 2 4C 3 



C4 : 



7-1 3P(T 3 



. C 2 +37 o /4 „ 2C7 4C > 
c 2 = -17^—, --t-w - -C 3 + =-, -rz + 4C + — — uj 2 



3Pct(7-1) V 3 P(7-l) 9(7-1) 
aC 2 7 20 2 16 

p(7^i) + y c<T+ y 

w 5 / 2C 4 , \ 4 

4CV +T77 i^ + ^ + i^V 3 - CB2) 



3(7-1) P(7-l) 3 

The dispersion relation is fourth order in the wave numb er fc, and fifth order in the angular frequency uj. It is 
identical to the dispersion relation given in equation 18 of I Johnson fc Kleinl ((2010); although they use co-moving 
radiation energy density and radiation flux for their analysis, there is no difference for the linearized e quations. The 
asymp totic behavior of the di s persio n relation for different parameters a and P are discussed in detail bv lBogdan et all 
(|1996|) and iJohnson fc Kleinl (|2010f ) and will not be repeated here. The analysis is done for a static background 
medium. Dispersion relation for a moving background medium is given bv lLowrie et al.l ()1999l ). In particular, they 
point out that the radiation hydrodynamic system is not Galilean invariant. The damping rate of some radiation 
dominated modes depends on the flow velocity, because the co-moving radiation flux seen by the gas depends on the 
flow velocity. 

In order t o check that the dispersion relation we get is consistent with previous results, we reproduce Figures 101.4 
and 101.5 of iMihalas fc Mihalasl (|1984l ) in Figure lB22l Those solutions are for boundary value problems, which means 
we assume a real number of an gular frequency uj and calculate a complex wave number fc. The dimensionless number 
Bo used in IMihalas fc Miha las (1984) is related to our dimensionless number P and C by equation ([25]) and another 
dimensionless number r is determined by the following equation 

^ ^ (B3> 

The Boltzmann number Bo is a dimensionless number which measures the importance of energy flux carried by the 
radiation field. When Bo — > 00, there is no energy exchange between radiation field and material. Equation [25] 
shows that the combination of our dimensionless numbers 1/(PC) has the same physical meaning. P, or equivalent 
r, is a dimensionless number which measures the importance of momentum exchange between the radiation field and 
material, as discussed in section RTTl 

The dispersion relation given above must be solved using numerical root finding methods. For convenience, we give 
four numerical examples of eigenmodes that span the parameter space defined by P and a which can be used by others 
to test their own codes In all cases, C = 10 4 and fc = 2-7T. 



DISPERSION RELATION FOR LINEAR WAVES IN RADIATION MHD 

In this appendix, we derive the dispersion relation for linear waves in radiation MHD used for the tests in 
section 16.11 As in Appendix [B] we assume a uniform, static background state, with a uniform magnetic field 
Bq = (B X fi, Byfi, B Zi q). We assume the the wave propagation direction is along x axis but the background 
magnetic field Bo can be along any direction. In general, v x , v y and v z are all non-zero. The analysis pro- 
ceeds as outlined in section [5] and will not be repeated here. In terms of the perturbed primitive variables 
SW — (Sp, Sv x , Sv v , Sv z , SP, SE r , SF r , x , 5F r ^, SF r , z , SB SB Z ) T , the linearized equations can be written as ASW = 0, 
where the coefficient matrix A is given in equation IC1I The dispersion relation for radiation MHD waves is given by 
det(A) — 0. This gives an eleventh order polynomial for u and a tenth order polynomial for fc. The solution can be 
simplified by assuming 2D solutions, in which Bo, the wave vector, and all the perturbed quantities are on the same 
plane. In this case, the dispersion relation will be reduced to an eighth order polynomial for uj and a seventh order 
polynomial for fc. Th e dispersion relation for radiation MHD waves including gravity in the short wavelength limit is 
discussed in detail bv Bl aes fc Socratesl (|2003l ). 

For convenience, we give four numerical examples of eigenmodes for both slow and fast magnetosonic modes that 
span the parameter space defined by P and a which can be used by others to test their own codes. In all cases, C = 10 4 
and fc = 2ir. 
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TABLE Bl 

Examples of solutions to the dispersion relation 1BH 



p 




0.01 






0.01 




a 




0.01 






10 




Sp 




10- 3 






10" 3 




Sv 


1.27177 x 10" 


3 +8.15409 x 10" 


5 i 


1.00012 x 10" 


3 + 6.91295 x 10" 


6 i 


SP 


1.61075 x HT 


3 +2.07402 x 10" 


A ■ 

i 


1.00019 x 10" 


3 + 1.36663 x 10" 


5 i 


SE r 


1.79137 x 10" 


8 +8.56498 x 10- 




6.64619 x 10" 


7 + 4.83820 x 10" 


5 . 


5F r 


-1.32035 x 10 


" 6 + 3.88814 x 10 


— R ■ 


—9.99976 x 10 


" 6 + 1.40748 x 10 


_7 . 

'% 


U! 


7.99077 + 0.512336i 




6.28393 + 


4.34354 x 10" 2 i 




P 




100 






100 




a 




0.01 






10 




Sp 




10" 3 






10 3 




Sv 


1.00000 x 10" 


3 +8.93099 x 10- 


s i 


9.99947 x 10" 


4 + 1.07703 x 10" 


& i 


SP 


1.00000 x 10" 


3 + 1.57240 x 10- 


7 i 


9.99998 x 10" 


4 + 1.60499 x 10" 


7 i 


SE r 


-3.55851 x 10- 


13 + 6.41394 x 10 


-10j 


—6 60096 x 10 


-9 4- fi 41 Sfi7 x 10 


~ 7 i 


SF r 


-1.00000 x 10" 


~ 9 + 2.10690 x 10" 


- 13 i 


-1.00130 x 10" 


- 9 + 5.35966 x 10" 


- u i 


w 


6.28319 + 


5.61151 x 10" 4 i 




6.28285 + 


6.76716 x 10" 2 i 




p 




100 






1 




(7 




0.1 






10 




Sp 




10" 3 






10" 3 




Sv 


1.00000 x 10" 


3 + 1.15555 x 10- 


7 i 


1.00000 x 10" 


3 + 3.32795 x 10" 


7 i 


SP 


1.00000 x 10" 


3 + 1.73114 x 10- 


8 i 


1.00000 x 10" 


3 + 2.94229 x 10" 


7 i 


SE r 


-1.0093 x 10- 


12 + 6.41394 x 10 


~ 9 i 


3.419000 x 10" 


10 + 1.11408 x 10 


-H 


SF r 


-1.00000 x 10" 


" 9 + 5.51807 x 10" 


- 13 i 


-1.00000 x 10" 


" 7 + 1.22263 x 10" 


- 10 i 


w 


6.28319 + 


7.26052 x 10" 4 i 




6.28319 + 


2.09101 x 10" 3 i 
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TABLE C2 

Examples of solutions for radiation modified slow mode 



p 
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4 i 


8 77920 x 10 -4 -1- 5 58924 x 10" 


6j 


sp 


1.52579 x 10" 


3 + 2.97056 x 10" 


4 i 
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